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, We consider the radiative feedback processes that operate during the formation of the first stars, 

■ including the photodissociation of H2, Lyman-a radiation pressure, formation and expansion of an H II 

region, and disk photoevaporation. These processes may inhibit continued accretion once the stellar mass 
has reached a critical value, and we evaluate this mass separately for each process. Photodissociation 
of H2 in the local dark matter minihalo occurs relatively early in the growth of the protostar, but we 
CO \ argue this does not affect subsequent accretion since by this time the depth of the potential is large 

enough for accretion to be mediated by atomic cooling. However, neighboring starless minihalos can 

1 (— 1 '' be affected. Ionization creates an H II region in the infalling envelope above and below the accretion 

disk. Lyman-a radiation pressure acting at the boundary of the H II region is effective at reversing infall 
from narrow polar directions when the star reaches ^ 20 — 30Mq, but cannot prevent infall from other 
directions. Expansion of the H II region beyond the gravitational escape radius for ionized gas occurs at 
masses ~ 50 — 100A/q, depending on the accretion rate and angular momentum of the inflow. However, 
^ ' again, accretion from the equatorial regions can continue since the neutral accretion disk has a finite 

thickness and shields a substantial fraction of the accretion envelope from direct ionizing fiux. At higher 
stellar masses, ~ 14OM0 in the fiducial case, the combination of declining accretion rates and increasing 
C^ ' photoevaporation-driven mass loss from the disk act to effectively halt the increase in the protostellar 

mass. We identify this process as the mechanism that terminates the growth of Population HI stars 
(i.e., stars with primordial composition) that have not been affected by prior star formation (Population 
III.l stars). We discuss the implications of our results for the initial mass function of these stars. In 
the Appendix we develop approximate solutions to a number of problems relevant to the formation of 
the first stars: the effect of Rayleigh scattering on line profiles in media of very large optical depth; the 
intensity of Lyman-a radiation in very opaque media; an approximate determination of the radiative 
acceleration in terms of the gradient of a modified radiation pressure; the determination of the flux of 
radiation in a shell with an arbitrary distribution of opacity; and the vertical structure of an accretion 
disk supported by gas pressure with constant opacity. 

Subject headings: stars: formation — early universe — cosmology: theory 



L INTRODUCTION 

There has been substantial recent progress in our theoretical understanding of how the first stars formed (Bromm & 
Larson 2004). In marked contrast to the case for contemporary star formation, the initial conditions for the formation of 
the first stars are believed to be relatively well understood: they are determined by the growth of small-scale gravitational 
instabilities from cosmological fluctuations in a cold dark matter universe. The flrst stars are expected to form at redshifts 
2; ~ 10 — 50 in dark matter halos of mass ~ lO^M© (Tegmark et al. 1997). In the absence of any elements heavier than 
helium (other than trace amounts of lithium) the chemistry and thermodynamics of the gas are very simple (Abel, Bryan, 
& Norman 2002, hereafter ABN; Bromm, Coppi, & Larson 2002, hereafter BCL02). There are no dust grains to couple the 
gas to radiation emitted by the protostar. There are no previous generations of stars to roil the gas out of which the stars 
form, nor is there any radiation other than the cosmic background radiation. Existing calculations have assumed that 
there were no significant primordial magnetic fields, thereby eliminating a major complication that occurs in contemporary 
star formation. However, even in the absence of significant primordial fields, it is possible that magnetic fields could have 
been generated in the accretion disk surrounding a primordial protostar (Tan & Blackman 2004), although even in this 
case the magnetic fields become dynamically important only after the star formation process is well underway, and they 
do not affect the initial conditions. Given this relative simplicity, there is some confidence in the results of numerical 
simulations that have followed the collapse of cosmological scale perturbations down to almost stellar dimensions (ABN; 
Bromm, Coppi, & Larson 1999; Yoshida et al. 2006; O'Shea & Norman 2007). This confidence is strengthened by the fact 
that it appears to be the microphysics of H2 cooling that determines the types of baryonic structures that are formed, 
and not, for example, the details of the initial power spectrum of fluctuations in dark matter density. The results of these 
simulations suggest that the initial gas cores out of which stars form are quite massive, Mcore ~ 100 — IOOOMq. 
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The observational imprint of the first stars depends on their mass: These stars were hkely to be of critical importance 
in reionizing the universe, in producing the first metals, and in creating the first stellar-mass black holes. The number 
of ionizing photons emitted per baryon depends on the stellar mass for m* < SOOMq (Bromm, Kudritzki, & Loeb 2001). 
The hardness of the radiation field is also sensitive to the mass (e.g. Tumlinson & Shull 2000; Schaerer 2002), so that 
He reionization can be affected. The effectiveness of the first stars in enriching the intergalactic medium (IGM) with 
metals and in producing the first stellar-mass black holes also depends sensitively on the mass of the star. A potential 
simplification in assessing these effects is that massive primordial stars are thought to have much smaller mass-loss rates 
than contemporary massive stars (Kudritzki 2002), so that the mass at core collapse should be quite similar to the initial 
mass. However, Meynet, Ekstrom, & Maeder (2006) have argued that if rotation is allowed for, then mass loss can be 
significant. Hegcr & Wooslcy (2002) showed that stars exploding as supcrnovae above about 260AfQ and between 40 and 
14OM0 should collapse directly to black holes, and they argued that such stars would providing relatively little metal 
enrichment. However, Ohkubo et al. (2007) followed the collapse and explosion of SOOM© and IOOOMq stars in two 
dimensions and concluded that substantial amount of metals could be ejected; they proposed that such supernovae could 
produce intermediate-mass black holes. Heger & Woosley (2002) also showed that for 140 < m* < 26OM0, the pair 
instability leads to explosive O and Si burning that completely disrupts the star, leaving no remnant and ejecting large 
quantities of heavy elements. Such supernovae produce a dramatic odd-even effect in the nuclei produced. Stars below 
~ 4OM0 are expected to form neutron stars, with more normal enrichment rates. In principle, metallicity determinations 
from high-redshift absorption line systems (Schaye et al. 2003; Norman, O'Shea, & Paschos 2004) or from very metal 
poor local stars (Beers & Christlieb 2005) can constrain the initial mass fimction (IMF) of the early generations of stars. 
Indeed, based on observations such as these, Daigne et al. (2004) argue that the stars responsible for reionizing the universe 
mostly likely had masses < IOOMq, and Tumlinson (2006) concludes that stars above 14OM0 could have contributed at 
most 10% of the iron observed in extremely metal poor stars (those with [Fe/H] < 10~^ — Beers & Christlieb 2005). 

In discussing the formation of the first stars, the terms "first stars" and "Population III stars" are often used inter- 
changeably, but this can lead to confusion. To be precise, we shall follow the conventions suggested by one of us at the 
First Stars HI conference (O'Shea et al. 2008) and define Population HI stars as those stars with a metallicity sufficiently 
low that it has no effect on either the formation or the evolution of the stars. The value of the critical metallicity for star 
formation — i.e., the value below which the metals do not influence star formation — is uncertain, with estimates ranging 
from ^ 10^^ Zq if the cooling is dominated by small dust grains that contain a significant fraction of the metals (Omukai 
et al. 2005) to ~ 10~^-^Zq if the cooling is dominated by the fine structure lines of C and O and there is negligible 
H2 (Bromm & Loeb 2003); if H2 cooling is included, Jappsen et al. (2007) argue that there is no critical metallicity for 
gas-phase metals. It is possible that values of the metallicity even less than 10~^Zq could significantly affect the evolution 
of primordial stars (Meynet, private communication). Among Population III stars, we distinguish between the first and 
second generations, termed Population III.l and III. 2, respectively: The initial conditions for the formation of Population 
III.l stars are determined solely by cosmological fluctuations whereas those for Population HI. 2 stars are signiflcantly 
affected by other stars. It is likely that Population III.l stars have a primordial composition, since it is hard to see how 
the gas out of which they form could be contaminated by even small trace amounts of metals without having been affected 
by radiation from the star that produced the metals. Stars affect the initial conditions for the formation of Population 
HI. 2 stars primarily by radiation, both ionizing radiation and Lyman- Werner band radiation that destroys H2 molecules. 
The latter reduces the cooling efficiency of the gas, allowing compression to heat up to the point that it begins to become 
coUisionally ionized; shocks associated with H II regions and supernova remnants can also ionize the gas. Once the gas 
has been partially ionized, HD cooling can become important, reducing the characteristic star-formation mass (Uchara & 
Inutsuka 2000). Greif & Bromm (2006) use the term "Population II. 5" to describe stars that form from gas in which HD 
cooling is important, whereas in our terminology these would be Population HI. 2 stars, but we believe that it is better 
to describe all essentially metal-free stars as "Population III." It should be noted that our definition of Population HI. 2 
stars includes all Population III stars that were significantly affected by previous generations of star formation, even if 
that did not result in significant HD production. Those Population III. 2 stars that form out of gas that is enriched in 
HD will typically be less massive than Population III.l stars — by about a factor 10 according to Greif & Bromm (2006). 
They infer that Population III.l stars are relatively rare, constituting about 10% by mass of all Population HI stars. 

In this paper we wish to estimate the characteristic mass of the first generation of stars (Population III.l). Even if they 
are relatively rare, they are critical in setting the initial conditions for the star formation that followed, and therefore in 
determining the reionization of the universe and the production of the first metals and the first stellar mass black holes. 
For contemporary star formation, it is believed that the IMF is set by a combination of gravitational fragmentation in 
a turbulent medium (Elmegreen 1997; Padoan & Nordlund 2002) and feedback effects. The characteristic stellar mass 
is of order the Bonnor-Ebert mass, m-BE T^/^/p^/^. However, not all of the initial core mass is incorporated into the 
final star, since contemporary protostars have powerful outfiows that eject some of the core mass (Nakano, Hasegawa, 
& Norman 1995; Matzner & McKee 2000). There are a number of feedback effects that occur in contemporary massive 
star formation (Larson & Starrfield 1971), particularly radiation pressure on dust and photoionization associated with 
the growth of an H II region. It remains unclear whether the upper limit on the contemporary IMF is set by feedback or 
by instabilities that afflict very massive stars. 

The clumps out of which the first stars form have total masses, including dark matter, of order 1O^M0; these objects 
are typically referred to as minihalos. Cooling by trace amounts of H2 generally leads to the formation of a gravitationally 
unstable core of baryonic mass 1O^~^M0 (ABN). In contrast to contemporary star-forming regions, the turbulence 
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in this gas is subsonic, due to the higher temperatures and the lack of internal and external sources of turbulence. As 
a result, gravitational fragmentation is much less effective — indeed, numerical simulations show no evidence for it (e.g., 
ABN, Yoshida et al. 2006, but see Clark, Glover, & Klessen 2008), and analytic calculations (Ripamonti & Abel 2004), 
including those that consider disk fragmentation (Tan & Blackman 2004), show no evidence for fragmentation either. It 
therefore appears that the mass of the first stars is likely to be set by feedback effects. 

Feedback effects can be classified as either radiative or kinetic. Kinetic feedback includes protostellar outflows and 
main-sequence winds. In contemporary star formation, outflows are believed to be hydromagnetically driven; in this 
paper, we assume that the magnetic fields associated with the protostar are too weak or too tangled to drive a strong 
outflow (see Tan & Blackman 2004 for a more extensive discussion of the effect of these outflows). Due to the absence 
of metals, main-sequence winds are very weak in the absence of rotation (Kudritzki 2002), although they could become 
important in the later stages of evolution if the star is rapidly rotating. Since we are primarily interested in the early 
stages of evolution of the star, we neglect kinetic feedback (i.e., outflows and winds). 

Radiative feedback includes radiation pressure, photoionization heating, and photodissociation. Radiation pressure can 
be due to continuum radiation or to resonance line scattering; the continuum radiation pressure can be due to electron 
scattering or to photoionization. Photoionization also leads to heating, which unbinds gas beyond the gravitational radius 
Tg = Gm*/Cg, where Cg is the isothermal sound speed of the gas. If the gas is initially in a disk, this process is termed 
photoevaporation. Finally, photodissociation destroys II2 the dominant coolant in neutral primordial gas. 

Most previous work has focused on the effects of continuum radiation pressure and photoionization heating in limiting 
the mass of primordial stars. Omukai & Palla (2001; 2003) focused on electron scattering, which leads to the Eddington 
limit on an accreting mass. For the case of spherical accretion at a rate of 4.4 x 10^'^ Mq yr^^, radiation pressure first 
becomes important at around 80Mq, leading to a dramatic swelling of the stellar surface. This, however, is a transient 
effect because an important part of the luminosity is due to accretion, which is reduced by the increase in the stellar 
radius. Only at masses around 300Mq does the internal luminosity become very close to the Eddington value, leading 
to runaway expansion of the star and, presumably, the end of accretion. Omukai & Palla (2003) considered a range 
of accretion rates. They found that if the accretion rate is smaller than 4.4 x 1O~^M0 yr^^, then the total luminosity 
remains sub-Eddington and the star continues to grow along the main sequence to arbitrarily large masses. On the other 
hand, if the accretion rate is somewhat larger than this critical rate, the Eddington limit becomes important at around 
1OOM0. Accretion at a rate based on the settling motions in the core of ABN is slow enough that the Eddington limit 
does not affect the final mass. However, these models ignored the influence of other protostellar feedback mechanisms on 
the infalling envelope. These models also assumed spherical symmetry, which leads to much larger photospheric radii and 
thus a softer radiation field than in the more realistic case of disk accretion (see §7 in Paper I). 

Omukai & Inutsuka (2002) considered the the combined effects of photoionization heating and continuum radiation 
pressure due to photoionization in the infalling envelope. They show that there is a critical stellar mass at which the 
hydrogen ionizing luminosity is sufficient to create an H II region, which rapidly spreads to large distances where its 
thermal pressure becomes dynamically important in slowing infall. However, the ionizing radiation force decelerates the 
inflowing gas, raising the gas density and therefore reducing the radius of the H II region. For spherical inflow, this 
mechanism is so effective that the radius of the H II region remains well below the gravitational radius r^, stopping any 
mass loss. They concluded that with this effect, there was no limit on protostellar masses below IOOOM0. Without this 
radiation force, they predicted a mass limit of order 300 A/q. As we shall see below, these conclusions are sensitive to the 
assumption of spherical accretion. 

In Paper I (Tan & McKee 2004) we modeled the growth of a primordial protostar from very small masses to large. We 
included the effects of rotation of the infalling gas, which led to the formation of an accretion disk around the protostar. 
The goal of this paper is to determine when the energy output from the protostar is sufficient to halt accretion and set 
the final stellar mass. This is an extremely complicated problem, the full solution of which requires three dimensional 
hydrodynamical simulations that include the generation, propagation, and dynamical influence of radiation. Furthermore 
these simulations must resolve a large range of scales: the protostar is of order 10 Rq, while the size of the quasi- hydrostatic 
core that encloses ^ lOOOM© is of order 1 pc, several million times greater. The demands on the time scale are even 
greater: the simulation may have to follow the evolution of the star over its lifetime '--^ 2 — 3 Myr (Schaerer 2002) while at 
the same time following the dynamics of an accretion disk with a characteristic time scale as short as 10^ s. The numerical 
simulations of ABN were able to resolve an even greater range of radii, but it will be some years before it is possible to 
meet the required dynamical range in time scales. As a result, we shall develop simple analytic models for the feedback 
interaction that we hope will provide a useful first step. 

We begin our discussion with a review the results of Paper I in 33 Feedback effects are then considered in the 
approximate order in which they become manifest. In fj3l we briefly discuss the effects of photodissociation. Lyman-a 
radiation pressure feedback, is considered in fj4l ^'^'^ in several Appendixes. The feedback from ionizing photons that 
can create an H II region is considered in fJSI After discussing shadowing by accretion disks in SJH and an Appendix, 
the feedback due to disk photoevaporation is considered in 33 Figure [1] gives an overview of these feedback processes 
occurring near the protostar. Finally, our conclusions are summarized in 33 

2. REVIEW OF PAPER I: PROPERTIES AND EVOLUTION OF PRIMORDIAL PROTOSTARS 

The radiative output from a protostar depends on the temperature and luminosity of its emitting components, which are 
the star itself (stellar photosphere), the boundary layer of the accretion disk with the star, and the larger scale accretion 
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Fig. 1. — Overview of the accretion geometry and feedback processes involved in primordial star formation, (a) Top left: Cross section 
of the accretion geometry: the dashed lines show streamlines of the rotating, infalling gas, with figure of revolution from each streamline 
separating 10% of the total infall from this hemisphere. The aspect ratio of the accretion disk is realistic, while the size of the star has been 
exaggerated for clarity, (b) Top right: The shaded region around the star shows the extent of the H II region, which at this relatively early 
stage is still confined inside the gravitational radius for the escape of ionized gas, Vg. Lyman-o radiation pressure feedback should be strong 
enough to prevent accretion in the polar directions, (c) Bottom left: The stellar mass and ionizing luminosity have grown, and the H II region 
is just in the process of breaking out of the accretion flow. Once a significant volume beyond rg is ionized, accretion from these directions is 
expected to be shut off. (d) Bottom right: Final stage of accretion involves shadowing of the equatorial region from stellar ionizing flux by 
the disk, which at the same time is photoevaporated. The competition between this photoevaporative outflow and the residual accretion rate 
sets the final mass of the star. 



disk. The luminosity of the star depends mostly on its mass. The size of the star then determines its temperature. The 
size of the star and the accretion rate determine the radiative properties of the both the boundary layer and the accretion 
disk, since emission from the latter is dominated by the inner regions. 

The size of the star depends on the accretion rate during its formation history. At lower masses there is a balance in 
the size set by the need to radiate the luminosity, which is mostly due to accretion, with a photospheric temperature that 
is likely to be close to ~ 6000 K because the opacity due to H~ rapidly declines in this temperature regime. Under the 
assumption of spherical accretion, Stabler, Palla, & Salpeter (1986) found the protostellar radius to be 

nc^90.8ml-fml% Rq (to,,2<0.1), (1) 

where m,^2 = m*/(100MQ) and m,^. -3 = w*d/(10~^ Mq yr~^). For the accretion rates typical of primordial star 
formation we see that the size is very large. For more massive protostars there is a transition once the star is about 
as old as its local Kelvin-Helmholtz time, then contraction proceeds towards the main sequence size, where accretion 
can continue. In Paper I, we found that for typical conditions, the protostar reached its main sequence radius at about 



2008.7.3.1900: DRAFT 



5 



1OOM0. According to Schaerer (2002), this radius is 

, 2 



r* ~ 4.3m° |^ Rq (main sequence) (2) 



to within 6% for 0.4 < m*,2 < 3. 

Thus almost ah the radiative steUar properties depend on just two parameters: the stellar mass and the accretion 
rate. Note that in principle these properties also depend on the angular momentum of the infalling gas, since if there 
was no rotation, then spherical accretion implies very high gas densities near the protostar, affecting the location of its 
photosphere. However, for any realistic amount of angular momentum, a disk forms whose size is much larger than r*, 
and the star's properties no longer depend on the rotation of the core from which the star forms. 

The accretion rate of Population III protostars depends on the density structure of the gas core at point when the star 
starts to form. This density structure is set by the balance of thermal pressure and self-gravity, which in turn depends 
primarily on the cooling properties of molecular hydrogen. This cooling creates almost isothermal cores at T ~ 200 K 
with an outer bounding density of about 10^ cm~^, which is the critical density of H2 cooling transitions (for H2 molecules 
interacting with atomic H). In fact the temperatures increase to several hundred K in the inner parts of the core because 
of the reduced cooling efhciency above the critical density. These basic features have been confirmed by numerical studies 
(ABN; BCL02; O'Shea & Norman 2007). The trigger for dynamical collapse is thought to be the rapid formation of II2 by 
three-body collisions at high densities ~ 10^" cm^'^, since this then dramatically increases the cooling rate in this region. 

ABN carried their calculations almost to the point of protostar formation, and at this time gas was flowing inward 
subsonically almost everywhere (except for O.lAf0 < M < IM©, where the inflow was slightly supersonic). Shu's (1977) 
expansion wave solutions for protostellar accretion are based on the assumption that the inflow velocity at this time is 
zero. Hunter (1977) generalized these solutions and showed that there is a discrete set of self-similar solutions that begin 
at rest att — —oo and have a constant infall velocity at the time of protostar formation {t ^ 0). One of these solutions, the 
Larson-Penston solution (Larson 1969; Penston 1969), has supersonic inflow (Mach number = 3.3 at t = 0; this solution 
is clearly inconsistent with the numerical results. In fact, the accretion flow appears to be a settling solution regulated 
by H2 cooling. Only one of Hunter's solutions corresponds to mildly subsonic inflow (Mach number =0.295 at t = 0), 
comparable to that found by ABN, and this is the solution adopted in Paper I. This solution has a density that is 1.189 
times greater than a singular isothermal sphere at t — 0, and the accretion rate is 2.6 times greater. 

Hunter's (1977) solution applies to an isothermal gas. Oniukai & Nishi (1998) and Ripamonti et al. (2002) have 
numerically calculated accretion rates for primordial protostars, and showed that the accreting gas is isentropic with an 
adiabatic index 7 ~ 1.1 due to H2 cooling; i.e., each mass element satisfies the relation P = Kp'' with the "entropy 
parameter" K = const. In hydrostatic equilibrium, such a gas settles into a polytropic configuration, which in general 
has P{r) — Kpp{r)'^^ . For an isentropic gas, we have Jp = J and Kp = K. In Paper I, we presented an analytic model for 
the protostellar accretion rate for isentropic gas. We allowed for the existence of an accretion disk around the protostar 
with a significant fraction of the stellar mass, 

= TO* + rrid = (1 + /d)m,, (3) 
with a fiducial value for the disk mass fraction fd = \. Following McKee & Tan (2002, 2003), we wrote the accretion rate 
as 

m^d^fP*— — , (4) 

where i^* is a numerical constant of order unity and iff = (37r/32Gp)^/^ is the free-fall time of the gas (measured at t = 0). 
For gas that is in hydrostatic equilibrium at t = 0, McKee & Tan (2002) showed that ^* ~ 1.62 — 0.96/(2 — 7^) to within 
about 1% for < 7p < 1; we have since verified that this is valid for 7^ < 1.2. To our knowledge. Hunter's self-similar 
solutions starting at t = —00 have not been generalized to the non-isothermal case. Q In Paper I, we therefore assumed 
that the accretion rate for the 7=1 case is 2.6 times greater than that for hydrostatic initial conditions, just as in the 
isothermal case. 

Feedback from the star, whether due to winds, photoionization, or radiation pressure, can reduce the accretion rate 
onto the star. We define a hypothetical star-disk mass, to^^^Oj and accretion rate, TO*d,Oi in the absence of feedback. In 
this case, the star-disk mass equals the mass of the core out of which it was formed, TO*d,o = M{r). The instantaneous 
and mean star formation efficiencies are 

£*d=-. , (5) 

m^d^o 

m^d m^d 

TO*d, M 

In our previous work, we assumed that the star formation efficiency was constant, so that e^d — ^*d- in the present work, 
we shall find that significant feedback does not set in until the star is fairly massive, so that we must distinguish the 

^ Fatuzzo, Adams, & Myers (2004) have given a comprehensive discussion of self-similar accretion solutions that start at t = 0, allowing for 
inflow velocity, overdensity relative to hydrostatic equilibrium, and non-isentropic gas (7 ■yp). Although they do not treat the time prior 
to protostar formation, their isothermal results for i > are consistent with Hunter's, as expected. For the non-isothermal case, Fatuzzo et 
al. (2004) present results for gas that is inflowing at r ^ 00, but these are not self-similar in that the accretion rate depends on where the 
integration begins (F. Adams, private communication). However, it is possible to generalize their treatment so that the Mach number of the 
inflow is constant. Presumably the overdensity and infall Mach number of the 7 = 1.1 analog of the mildly subsonic Hunter solution axe 
similar to those of the isothermal solution; if they are exactly the same, then the accretion rate would be about 2.0 times that for the case of 
hydrostatic initial conditions, somewhat less than the isothermal value of 2.6. 
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instantaneous and mean values. In Paper I, we found that the accretion rate onto the star-disk system is 



m^d = 0.026 

where 



{m/MqYI\ 



Mq yr-\ (7) 



~ 1.88 X 1012 cgs VSOOK/V «H / ^ 

is a measure of the entropy of the accreting gas. Here T^g = T + M'^'^urb/'^ effective temperature that includes 

the effect of turbulent motions; we have added a prime to the Toff defined in Paper I to distinguish it from the effective 
temperature of a radiating atmosphere. Expressing the accretion rate in terms of the stellar mass, which equations p|) 
and ([6]) imply is = e*c;Af/(l + fd), we find 

e^d^K'^'/' 



m^d = 0.026 



(l + /d)3A(m,/Mo)3A 



Mq yr-\ (9) 



With K' = e^,d — e*d = 1, this result is in good agreement with the results of Omukai & Nishi (1998) and Ripamonti et al. 
(2002); since their ID calculations did not allow for disks, this comparison is made for fd — 0. Note that this agreement 
validates our use of the Hunter mildly subsonic solution to infer the accretion rate. If e^d ~ 1 (i-e., no feedback) and 
fd — \ then the accretion rate onto the star + disk is 

/^/15/7\ 

m,rf,_3^3.20^7^ , (10) 

where henceforth it will be understood that numerical evaluations denoted by have = 1 and fd = \- In this case 
the accretion rate onto the star (which may be primarily through the disk) is 3/4 of this [since m* = m*d/(l + fd)]- 

Our estimate of the accretion rate is somewhat above that estimated by ABN, but this is to be expected since their 
calculation stopped prior to the formation of the protostar. Indeed, at the time at which the protostar first forms (i = 0), 
the accretion rate at any finite radius r [i.e. m(r) = 47rr2p|wr|] in a self-similar, isothermal collapse is smaller than the 
value it has at times t > r/cg, where Cg is the isothermal sound speed. Equivalently, at a given time, the accretion rate 
at radii r > Cgt is less than that at small radii, r ^ Cgt. In Shu's (1977) solution for the collapse of a singular isothermal 
sphere, the accretion rate at a given time is zero outside the expansion wave at r = Cgi; inside the expansion wave, the 
accretion rate smoothly increases to 0.975Cg/G as r — > 0. For the Larson-Penston solution, the accretion rate at a given 
time i > increases from 29Cg/G at large radii (r 3> Cgt) to 47c^/G at small radii (r ^ Cgt). For Hunter's mildly subsonic 
solution, which we have suggested is closest to the numerical simulations, the accretion rate increases from O.7OC5/G at 
large radii to 2.58c^/G at small radii (Hunter 1977), an increase of a factor 3.7. This demonstrates that caution should 
be exercised in inferring accretion rates at late times from those measured at early times, which is commonly done in 
simulations (e.g., ABN, Yoshida et al. 2006, O'Shea and Norman 2007). 

The age of the star when it reaches a mass m* is (Paper I) 

/I _L f \ 1"/''' / \ 1*'/''' 

ty^. = 27(iii^) K'-^^niai^] ^ 2.92 X lO^if'-is/ , (H) 



where ty-^ = t/{l yr) and where it is the mean star formation efficiency e^,d that enters. The resulting stellar mass is 

^ OmbK'^ h'^J Mq. (12) 

Bromm & Loeb (2004) have carried out a 3D simulation of the accretion onto the protostar for the first 10^ yr, and for 
K' ^ 1 our result is within a factor ~ 2 of theirs for this time period. (However, it should be noted that an extrapolation 
of their result to times beyond 5 x 10^ yr gives a mass less than half our estimate of the mass of the star plus disk. It 
remains to be determined whether such an extrapolation is valid.) 

With this estimate of the protostellar mass, it is possible to calculate the maximum possible mass a primordial star could 
have. Schaerer (2002) has calculated the main sequence lifetimes of primordial stars with no mass-loss for < 5OOM0, 
and his results are accurately described by the expression t^s — 2.7m~ j^** Myr for lOOM© < m* < SOOM©. If we assume 
that the accretion is not limited by any feedback processes (e*d = 1), that Schaerer's results can be extrapolated to higher 
masses, and that accretion during the relatively short post main-sequence phase is negligible, then we find 

/ 1 \0-86 

dt ~ 1900 K'^-^^ Mq ISOOif'^-^^ Mq. (13) 

V 1 + Jd / 

The maximum possible stellar mass is thus controlled by the value of the entropy parameter of the core. 

In Paper I we also considered the effect of rotation. Rotation of the infalling gas has a dramatic effect on the evolution 
of the protostar, since it leads to much smaller photospheric radii and correspondingly higher temperatures and accretion 
luminosities. We parameterized the rotation in terms of 

f — ^^rot(?'sp) 

/Kop = 7 T, (14) 

WKop(?'sp) 
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the ratio of the rotational velocity to the Keplerian velocity measured at the sonic point at rgp. ABN found /xop '^-^ 0.5 
independent of radius, so we adopt this as a fiducial value. We then showed that the accreting gas formed a disk with an 
outer radius 

cm, (15) 
cm. (16) 

3. PHOTODISSOCIATION FEEDBACK 

As the protostar grows in mass, it begins to emit copious amounts of non-ionizing ultraviolet radiation (far-ultraviolet, 
or FUV, radiation) , as shown in Figure [21 This radiation can photodissociate the H2 that is critical for cooling the 
accreting gas (Omukai & Nishi 1999; Glover & Brand 2001); its dynamical effects are considered in the next section. 
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Fig. 2. — Evolution of Lyman- Werner photon luminosity from the fiducial model of primordial star formation, including effects of stellar 
feedback. The total (solid line) and contributions from the protostellar surface (long-dashed line), boundary layer (dashed line), and accretion 
disk from r < lOr* (dotted line) are shown. 

Once the molecular coolants in the accreting gas are destroyed, the adiabatic index rises from 7~l.lt0 7= |. If the 
gas were able to continue contraction, it would heat up until the temperature became high enough {T ~ 10'* K) to excite 
the Ly-a transition. For T > 10^ K, the adiabatic index would then drop back to about 1. 

Can the protostar continue to accrete when 7 ~ |? If one considers the related problem of the gravitational stability 
of polytropic gas spheres, one might be led to conclude that accretion would stop: polytropic stars are stable against 
gravitational collapse for 7 > | (Chandrasekhar 1939), and the same is true for polytropic gas clouds even if 7p < | 
(McKee & Holliman 1999). However, there is a crucial distinction between collapse onto a protostar and the contraction 
of a gas cloud prior to protostar formation, and that is the presence of the central protostar, which is effectively a mass 
singularity. The stability analyses cited above assumed that there was no mass singularity at the origin. When one is 
present, the problem is analogous to that of Bondi accretion, the accretion of non-self gravitating gas onto a star; this can 
occur for 7 = f (e.g., Shapiro & Teukolsky 1983). The problem of protostellar accretion, which includes the self-gravity of 
the gas, has been considered by Fatuzzo et al. (2004) for a wide range of conditions. They showed that gravity dominates 
over pressure close to the protostar, so that accretion can occur, provided 7 < |. It is straightforward to see why: In 
supersonic inflow, the density scales as p oc r"'^/^, so that the temperature T oc r~^(''~^y^ rises more slowly than the 
kinetic energy per unit mass cx provided 7 < |. They demonstrated that the accretion rate for a singular, initially 
isothermal sphere with 7 = 1.6 is only slightly smaller than for the case in which 7=1. 

The argument of Fatuzzo et al. applies to the inner, supersonic region of infall. What about the outer, subsonic region? 
There the density varies as a higher power of radius (e.g., for Umfaii = const., p cx r~^), so that pressure can overcome 
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gravity at a lower value of 7. To see this more quantitively, consider the case of a singular isothermal sphere with 7 7^ 1. 
Assume that the initial density of the sphere is A times greater than the equilibrium value. If the gas is flowing inward 
at a velocity —v^o far from the protostar (i.e., at large values of the similarity variable x, which is just r/[cst\ in the 
isothermal case), then v varies as 

2(A-1) (4-37K0 
V = -Woo 5 1 ■ (17) 

X X'^ 

(Fatuzzo et al. 2004; we have corrected a typo in the last factor). Thus, for 7 > |, pressure forces will tend to decelerate 
the flow; however, this can be overcome by a suitable overdensity A. We have confirmed this by numerical integration of 
the equations given by Fatuzzo et al.: For Voo > and 7 < |, accretion is possible for A > 1; for | > 7 > |, accretion is 
possible provided A exceeds some threshold. For primordial star formation, we estimate 7p = 1.1 and Voo/cg — 0.3 — 0.5; 
accretion can occur in this case for 7 = | for A > 1.16, 1.28, respectively. These overdensities are quite modest (for 
example. Hunter's 1977 subsonic infall solution has Voo/cs = 0.295 and A = 1.189), so we anticipate that photodissociation 
should not prevent protostellar accretion. The value of the overdensity is likely to vary from one protostar to another, 
however, so it is possible that in some cases it would be too small to permit accretion. In such cases the infalling gas 
would decelerate; once it is stationary, however, it could resume accretion when it is overtaken by an expansion wave, as 
shown by Fatuzzo et al. Our numerical calculations show that the increase in 7 from 1.1 to | has only a minor effect 
on the accretion rate, diminishing it by less than 20%. We conclude that photodissociation of molecular coolants by the 
protostar does not have a significant effect on its accretion rate. 

On the other hand, collapsing cores that do not contain a protostar, but that are close enough to a protostar that their 
molecular coolants are destroyed, will cease collapsing if their central temperature is low enough (< lO'' K) that 7 exceeds 
|. Thus, FUV emission from the first stars is potentially quite effective at suppressing star formation in their vicinity. 
We can estimate the distance over which star formation is suppressed from the work of Glover & Brand (2001). As in 
Paper I, we assume that the core is in approximate hydrostatic equilibrium and is characterized by an entropy parameter 
K. We find that the time to dissociate II2 is less than the free-fall time tg if the core is within a distance 
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of the protostar, where 5lw is the photon luminosity in the Lyman- Werner bands, X2 is the fractional abundance of 
II2, /abs is the fraction of the Lyman- Werner flux absorbed by the II2, /diss is the fraction of absorptions that result 
in dissociation, and n is the mean density of H nuclei. Thus, even a lOOM© star, which has S'lw — 10^^ s~^, can 
suppress star formation in an existing core only if the core is relatively nearby. A more detailed analysis by Susa (2007) 
comes to the same conclusion. Ahn & Shapiro (2007) model both dissociation and ionizing feedback and also find a 
relatively weak suppression of PopIII.2 star formation by neighboring PopIII.l stars. Whalen et al. (2008) have presented 
multi-dimensional numerical simulations of these processes. 

4. LYMAN-a RADIATION PRESSURE 

The second feedback effect of FUV radiation is radiation pressure acting on the Lyman absorption series in the infalling 
neutral gas. This effect has been considered previously by Oh & Haiman (2002), who studied feedback effects in halos 
with virial temperatures above 10^ K, which are more massive than those we consider. They concluded that Lyman-a 
radiation pressure could be important, but did not find any constraint on the mass of the star that could form. Our 
work improves upon theirs in several respects: we include stellar continuum photons injected away from line center as 
well as Lyman-a photons emitted in the H II region; we allow for Rayleigh scattering; we include the limitations on the 
radiation pressure set by two-photon emission and by the blackbody constraint; and we allow for the effects of rotation 
in the infalling gas. 

Since conditions are very opaque, the Lyman-a radiation can be considered to be isotropic. The Lyman-a radiation 
pressure is then 

P. = = (19) 

where Ua and Jq are the energy density and mean intensity of the Lyman-a radiation. The estimation of Ja is complicated 
by the fact that Lyman-a photons can diffuse in frequency as well as in space, and that at the optical depths we are 
considering, the transfer is dominated by the damping wings of the line profile (Adams 1972). This problem is far too 
difficult to treat analytically for the geometry and dynamics that we are using to model the protostellar accretion. We 
therefore make the following substantial approximations when evaluating Ja at the outer boundary of the H II region, 
fmi (which may be at the surface of the protostar), and at a particular polar angle. (1) The axially symmetric geometry 
can be replaced by an equivalent slab geometry. The effects of spherical divergence are incorporated by normalizing the 
mean intensity to the flux at rnii- The slab column is set equal to 20% of that in the infalling gas based on the discussion 
in Appendix [Cl (2) The anisotropy in the optical depth can be accounted for by taking the harmonic mean of the opacity, 

1 1 f dA 



Teff A J T(r) 

(see Appendix [D|) . In practice, the escape of photons is primarily controlled by the minimum optical depth, which is in 
the polar direction, so in our numerical models we evaluate feff with a column density that is 20% of the column in the 
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vertical direction from the point of interest. (3) Finally, we assume that the effect of the velocity field can be approximated 
by a Doppler line profile of suitable width (see below) . 

The propagation of resonance line photons in very opaque media has been treated by a number of authors (Adams 
1972; Harrington 1973; Hummer & Kunasz 1980; Ncufeld 1990). Let the mean optical depth in the line be 

1 



ydAv. (21) 



In terms of the normalized frequency x = Av/Ai'o, we have — — f(j)x^ where <f>x is the line profile. In the Doppler 
core, the line profile is ~ exp(— a;^)/-y/7r; in the damping wings it is — where a is the ratio of the natural 

line width to the Doppler width. In applications, we generally have a ^ 1, and in that case the optical depth at line 
center, tq, is related to mean optical depth by tq = t(/)o = t/^tt. The opacity is = k4>x, and the mean free path is 

— l/^a;. 

As shown by Adams (1972), resonance photons escape in a single longest excursion from line center. After n scatterings, 
the escaping photon has a frequency shift Xe — n^^'^ and it has traveled a distance rt^/^^e — n^^'^/ (Kcjje), where £e — i{xe), 
etc. In order for the photon to escape, this distance must be the size of the region, L — tl/R. This implies n^/^ ~ tl^c 
and cCe — Tz,</>e, which in turn gives the characteristic frequency of the escaping photons as Xe (afx)^/'^. The total path 
length traversed by the escaping photons is about n^^^L. As a result, we expect the mean intensity to exceed the incident 
intensity by a factor of about n^'^ ^ (afi)^/"^. 

The velocity field has contributions from thermal motions, turbulent motions, and the overall flow. Thermal and 
microturbulent velocities are naturally included in Avd. In the cases we shall consider, the overall flow is highly opaque, 
so it generally does not contribute to the random walk of the photons. If the velocity width of the flow Avf (including 
macroturbulence) is small compared to the line width of the escaping photons, Av^ — (afcsy^^Avu (from Neufeld 1990), 
then the flow has negligible effect on the escape of the photons. On the other hand, if Avf ^ Ave, then the effective 
column density of the gas is reduced. For example, in the simple case in which the velocity varies linearly with the column 
density, photons will interact with only a fraction Ave/ Avf of the gas. If Aveo is the line width in the absence of any 
flow velocity, then one can show that the effective column density is reduced by a factor of about {Avgo/ Avf)^^^ from 
the total value. In our numerical models we always set Avd = 12.9km s~^, the value appropriate for T = 10'* K gas, the 
assumed temperature of the infalling neutral gas near the protostar. We set Avf equal to half the difference in radial 
velocities of the inner and outer edges of the slab. If Aveo > ^Vf, which is not usually the case, we reduce the effective 
column by the factor (Aueo/Au/)'^/^. 

Appendix |B] describes the general enhancement in the intensity of photons that are trapped by the Lyman-a damping 
wings, and, if the columns are very large, by the opacity due to Rayleigh scattering. Thus, photons from the protostellar 
continuum, outside the frequency interval defined by Xe, can contribute to the radiation pressure. The enhancement in 
intensity leads to an increase in the radiation pressure so that the momentum transferred from the radiation to the gas 
is F/c in each optical depth. As shown in Appendix iBl the isotropic component of the radiation pressure is 

p 4nJ,,o An f f 8.25N'J%Av-^%'iF^Jc) 1 

frad,iso = — 5 — ^ ^ I dvi Mm < By., — -jz ) , (22) 

3 3 7 \ Min[l, 2.&2Nll,%Av-^f] + 5.41[£f //(^,)] + m) J 

where A'cfr. 20 = A'cff/(10^° cm~^) and B,j. is the intensity of a blackbody with a temperature equal to that at the 
protostellar surface. When this limit is evaluated for the reprocessed Lyman-a photons, the intensity is limited to that of 
a blackbody at the temperature of the ionized gas in the H II region (see Appendix lB.3p . This expression is valid provided 
Toff > 1/a, corresponding to A'off > IQ^^Avjj g cm~^ for Lyman a. 

What is the condition for the radiation pressure to halt the accretion? We assume that the accreting gas is inside the 
sonic point, so that the gas pressure is negligible. For steady, radial flow, the equation of motion of the gas is 

^'^ _ ^-Prad, iso pGm^ 
dr dr ' 

since we have shown in Appendix [C] that the radiative force can be represented by the gradient of the isotropic component 
of the radiation pressure. We assume that the radiation pressure builds up rapidly over a distance small compared to the 
radius r; this is justified below. Then constancy of the mass flux implies pv ~ const. If the radiative force is to stop the 
flow in a small distance, then the gravitational force must be negligible in comparison. We then have 

^ (pf'+Prad,iso) ^0, (24) 

so that pv^ + Prad, iso — const in the deceleration region. When gas enters this region, the radiation pressure is small 
and u ~ tiff, the free-fall velocity; as the gas decelerates, the radiation pressure rises and v drops. The inflow will be 
halted if the radiation pressure at the inner edge of the deceleration region is Prad.iso = Ps'L's^ where pff is the density 
in the freely falling gasQ We have verified this simple argument by solving for the flow in the case that the flux varies 

^ Jijina & Adams (1996) have given an alternative criterion based on treating the radiative force per unit mass as the gradient of a potential. 
Their approach is appropriate when one knows the spatial variation of the force in advance, which is not the case here. One can show that the 
two approaches are equivalent if the radiative force per unit mass falls off rapidly with r. 
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as r"'^*', with kp > 2; such a faster than spherical fall-off in the flux is expected when the density distribution is not 
spherically symmetric, so that flux will escape into regions of lower opacity (as in the case of the "flashlight effect" — Yorke 
and Bodenheimer 1999). We find that the radiation pressure required to halt the infall is within about 10% of pevg for 
kp > 2.5. In order to reverse the inflow and eject the matter, the radiation pressure must be twice this, Prad.iso — 2pw^. 
Of course, a steady radial flow cannot reverse direction. To see that the factor 2 is required, one can imagine that the 
flow is inward over half the sky and outward over the other half; to maintain the same accretion rate, the density would 
have to be twice as large. If the gas is initially in a disk, there is no infall to start with and the pressure required for 
ejection is pv^ = pvs/2. 
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Fig. 3. — Protostellar mass scale at which Lyman-a radiation pressure becomes twice the ram pressure of the infalling gas at the edge of 
the H II region along the polar direction as a function of /kcp ■ At this point the radiation pressure is expected to reverse the accretion in the 
polar direction and evacuate a polar cavity, through which Lyman-a photons can escape. Thus this feedback mechanism will act to reduce 
the efficiency of accretion, but will not significantly impede the growth of the star, since most mass is accreted from directions away from the 
polar axis. 

We evaluate this criterion along the polar axis at the edge of the H II region, which is where the ram pressure of infalling 
gas is a minimum and where breakout should occur first (Figure ^ . At low values of /kop the breakout does not occur 
until the star has reached several hundred solar masses as the photosphere is very large and cool, producing little FUV 
flux. At reasonable values of /xep 0.1, polar breakout can occur relatively early, at ~ 20 Mq. This is the point in the 
protostellar evolution when the star is starting its rapid contraction to the main sequence, and the surface temperature 
and luminosity are thus rising. For these values of /kop the densities and ram pressures become significantly greater as the 
sight line moves away from the pole. By the time that the radiative flux from the star is large enough to reverse the flow 
in these directions, a polar cavity would have been blown out, thus reducing the enhancement in the radiation pressure 
due to trapping of photons. Thus although Lyman-a radiation pressure can act to reduce the efficiency of accretion, we 
expect it to be unable to stop it. Even the reduction of the accretion efficiency is likely to be relatively modest, since 
even a small polar cavity can dramatically reduce the radiation pressure in the H II region. In the following sections we 
consider other feedback mechanisms that are more effective at limiting accretion, although they do so at higher masses. 
When following the stellar evolution to these masses we shall assume the reduction in accretion efficiency due to Lyman-a 
feedback is negligible. 



5. IONIZING FEEDBACK AND BREAKOUT OF THE H II REGION 

5.1. Photoionization Heating 

Extreme ultraviolet {hv > 13.6 eV) radiation from the protostar can ionize infalling neutral gas, creating an H II 
region. The temperature of the ionized gas is about ~ 2.5 x 10** K, based on the models of Giroux & Shapiro (1996) 
and Shapiro, Iliev and Raga (2004) with stellar spectra. The thermal pressure P = pc^ of the ionized region is typically 
much greater than that in neutral gas of the same density because of the elevated temperatures and sound speeds: 
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a = 11.6(r,/104K)i/2kni The pressure gradients created at this ionized-neutral boundary can become steep enough 
to cause the H II region to expand and to dramatically reduce the accretion of gas to the star. In this section we attempt 
to calculate the point in the protostellar evolution at which this transition occurs. This problem has been considered 
previously by Omukai & Inutsuka (2002). The new feature in our treatment is that wc allow for the rotation of the 
infalling gas, which can significantly reduce the density near the protostar. As we shall see, this completely changes the 
evolution of the H II region. 

As in Paper I, we approximate the density distribution of the infalling gas by the Ulrich (1976) solution. The gas is 
assumed to be spherically symmetric far from the protostar, and each mass element conserves its angular momentum 
as it falls ballistically toward the star. Terebey, Shu, & Cassen (1984) showed that this solution matches on to an 
expansion-wave solution for the gravitational collapse of a singular isothermal sphere. The resulting density distribution 
is 

^ 4^r3/2(2Gm)i/2' ^^^^ 



where fi = cos 0, 

( 



1/2 



1 



JJ_ jj_ 

\ A^o / ^J.Q 

and iio is the value of far from the protostar. The two angles are related by 



^ (26) 



(27) 



r /io(l-Aio)' 

which shows that /iq > M- t^i^ converges toward the disk plane. Ulrich assumed that the disk had negligible mass, so 
that m = TO* in equation (|25p . In our case, m varies smoothly from m^d to as r shrinks from being much greater than 
Vd to being much less than rd. This variation in the mass acting on the infalling gas leads to small, unknown deviations 
from the Ulrich solution. In view of the necessarily approximate nature of the solution and the fact that m*^ and m* 
differ by only a factor | in the fiducial case, we shall set m = m^:d in applying equation (|25p . 

The density factor ip depends on both the current direction cosine, ji, as well as the initial one, /io, with the two 
being related by the cubic equation (j27[) . In our analysis, it is convenient to have an approximation for -0 in which the 
dependence on /^o is eliminated: 

1/2 ^ 

(28) 



1 + Max(^2/3^1_^) 



0.5 (C - 1 + 3|C - 1|) + 3Ai2/3Min(l, C) ' 

where C = r^/r. This is exact for all r ai 6 = where /i = /j,o = 1, and at 6* = 7r/2, in the plane of the disk. It is also 
exact aX r ^ Td for all 9. For r < it is accurate to better than 20%; for r > r^, it is accurate to within a factor 2. To 
simplify our results, we approximate this further for r < 0.5rd and take 

1/2 ^ 



This approximation is quite accurate at /i = | (better than 20% for r < rd), but it deteriorates away from there, 
underestimating the density by a factor ~ 2 in the plane for r = ^rd (the accuracy improves as r decreases). Nonetheless, 
it suffices to give an analytic estimate for the behavior of the H II region. 

5.1.1. Early Evolution of the H 11 Region 

H II regions are bounded by ionization fronts. Ionization fronts that move faster than about 2ci with respect to the 
neutral gas, where Ci is the isothermal sound speed of the ionized gas, are termed "R-type;" such fronts have little 
dynamical effect (Spitzer 1978). However, if the velocity of the front slows below 2ci, a shock forms in front of the 
ionization front and the velocity of the front into the shocked gas falls to ~ c^/2ci, where c„ <C Ci is the isothermal sound 
speed of the shocked neutral gas; such ionization fronts are termed "D-type." When the H II region first forms, it is 
embedded in gas falling inward with a velocity vg 3> 2ci. As a result, the ionization front is initially R-type, and the 
radius of the H II region, thii, is determined by ionization balance. 

Since the density of the infalling gas depends on the angle 9 relative to the axis of rotation, thii depends on angle also. 
We determine this radius in the sector approximation, in which ionizations balance recombinations in each element of 
solid angle: 

— — r'^a'^'^^UeUpdr, (30) 

ail in jj.^ 

where S is the rate of emission of ionizing photons and a'-^-' ~ 2.6 x 10~^3y^~0-8 (,]^3 g-i jg -j-j^g recombination rate to 
the excited states of ionized hydrogen. In writing equation pop . we have made three approximations. First, we have 
assumed that the rate of emission of ionizing photons is much greater than the rate of accretion of hydrogen atoms so 
that ionizations and recombinations are very nearly in balance (note that advection of neutral H into the H II region is 
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allowed for in the numerical models). For a mass accretion rate of 10"'^ Mq yr~^, the hydrogen accretion rate is about 
3 X 10^^ s^^. The mass at which the ionizing photon luminosity exceeds this value depends on /xep; for the fiducial case 
of /kcp = 0.5, this occurs at about to* ~ 20Mq. Second, we have assumed that in the outer parts of the H II region, 
where the helium is singly ionized for stellar temperatures ~ 10^ K, each recombination of He'^ results in one H ionization, 
whereas it actually results in about 2/3 of an ionization at the relevant densities (> 10'* cm^'^; Osterbrock 1989). In fact, 
the abundance of He is sufficiently small (~ 0.08) that we shall henceforth neglect it in our analytic estimates (however, 
we do not neglect its contribution to the mass density, nor is it neglected in the numerical calculations). Finally, we have 
ignored photoionization from the n — 2 level of H, so that our calculation somewhat underestimates the size of the ionized 
region, although this is not very important at the densities resulting from realistic values of /xcp- 
With the density distribution given by equation equation (|5D|) becomes 

S^f^.S^J, (31) 
87r^^GTO,d 

where — 2.20 x 10^^"* g is the mass per hydrogen and 

I^r^^dr. (32) 



We have set to = m^,d in equation (|25p in accord with the discussion following equation ((27|) . Equation ((3T|) reduces to 
that of Omukai & Inutsuka (2002) for -0 = ln(rHii/r'*) (and if is replaced by nip, m^^d by m* and m,d by to*). As 
shown by Omukai & Inutsuka, an H II region in an r"'^/^ density profile is confined to the vicinity of the star for S < S'cr 
and expands to exponentially large distances as S increases above S'cr. Numerically, we have 

/OCX 0.8 -2 /„c\0.8 7^/30/7 

S„ = 3.07xl05°(±^) ph ^ 2.36 X 10^1 (M) ^i— ph s-^. (33) 

By comparison, the ionizing luminosity of a Pop III star is 

5 ~ 7.9 X lO''^ (/)sTOi | Phs-i, (34) 

which for 05 = 1 is a fit to Schaerer's (2002) results for the ionizing luminosity of main sequence primordial stars; the fit 
is accurate to within about 5% for 60Mq < to* < 300Mq. As shown in Paper I, the ionizing luminosity can be less than 
the main sequence value (0s < 1) when the star is still contracting toward the main sequence, and it can be greater when 
it is accreting while on the main sequence; for the case illustrated in Paper I, 0s ^ 2. If the accretion rate is not reduced 
by feedback effects, S would not exceed S'cr until to* > 21hK'^^l^'^ Mq for T4 — 2.5. However, as we shall see, rotation 
makes the factor / small and allows the H II region to expand until it is almost as large as the disk even when the mass 
is less than this. 

At early times we have rnii ^ Td so that we can use approximation ()29|1 for the density. As a result, we find 

- / ~ 1 ( mil \ ^ 
Scr " ~ ^{\^ii)\rd ) ■ 
With equations ([33)l and (pi)) , we then obtain 

/ rp \ _1.25 

^ = 1.01(1 + + fd)'/Vj' (^) (36) 

Td V2-5/ TO*d,_3 



(35) 



, ^ X 0.4 ™47/28 

O.37(l+M)^/^0f (g 9^; (37) 



recall that indicates that we have taken to*^ — |to*. We see that for to*, 2 ^ 1 we have thii < 0-5rd, so that our 
approximation for the density, equation (|29p . is reasonably accurate at early times. The radius of the H II region is then 

.H„ = 5.4OxlO-(l + M)^/^0f cm (TO*,2<1), (38) 

where we have set the star formation efficiencies e*d and e^d equal to unity and where we have made the approximation 
83/28 ~ 3 in the exponent of to*, 2. As the stellar mass increases above IOOMq, the approximation for the density, 
equation (|29p . increasingly overestimates the density except near the equator; as a result the radius of the H II region, 
r-Hii , becomes larger than the value given in equation ([55]) except near the plane of the disk, where the high density traps 
the H II region. As remarked above, for S > S'cr, which occurs for to*, 2 > 2.75 if the accretion continues unabated by 
feedback, rnii increases exponentially with S (Omukai & Inutsuka 2002). 

5.1.2. Later Evolution of the H 11 Region: Suppression of Accretion 

According to equation (I38p . the radius of the H II region expands on the protostellar evolution timescale ~ 10* yr, 
which is far longer than the dynamical time rHii/2ci ^ 10^ yr. As a result, the velocity of the ionization front relative to 
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the infaUing gas is very nearly equal to the free-fall velocity. The first phase of evolution of the H II region ends when it 
expands to the point that the radius becomes comparable to the gravitational radius, 

rg = 2 = 3.92 X lO'^'^Edd 77^ m*d,2 cm, (39) 



T4 

where we have taken the gravitating mass to be m^d- Here we have allowed for the decrease in the radiation pressure due 
to electron scattering through the factor 

</.Edd = 1 - 7 , (40) 



;dd 



where Lsdd = 47rGmc/KThomson is the Eddington limit. In Paper I we found that L/LEdd ^ 0.6—0.8 for m — to* ~ 10^ Mq, 
which corresponds to 0Edd — 0.2 — 0.4. Equations (j38|l and ([39|) relate the protostellar mass to rmi/rg, 
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(41) 



_(l + Ai)l/4(l + /rf)lV280V4_ 

Keep in mind that this relation is valid only for m*^2 ^ 1, so that /kop cannot be small. This condition is satisfied insofar 
as the simulations of ABN are representative of the angular momentum of the accreting gas, since they give /kcp ~ 0.5. 

When the H II region expands to the point that Vff — 2ci, a shock front forms and the ionization front becomes D-type; 
this occurs at rnii = rg/2. The accretion rate through the H II region will begin to decrease at this point. Since the 
shocked neutral gas is denser than the ionized gas in the H II region, the accelerating expansion of the shocked shell is 
subject to the Rayleigh- Taylor instability, and as a result it is difficult to estimate by how much the accretion is reduced. 
While the shell is moving slowly, it can fall onto the disk and accrete that way. However, once the shocked shell is moving 
faster than the local free-fall velocity it seems unlikely than any significant further accretion can occur. To obtain an 
approximate upper limit on the accretion through the H II region, we assume that, from a given direction, the accretion 
is unimpeded until the H II region has expanded to a radius equal to rg . Because of the declining density distribution in 
the infall envelope, the H II region is expanding relatively rapidly at this stage and so soon ionizes a large region beyond 
Tg, which we expect leads to a substantial reduction in accretion rate from the affected directions. This approximation 
needs to be checked with multi-dimensional radiation- hydrodynamical simulations. It is important to bear in mind that 
this suppression of the accretion occurs only in the ionized gas. Gas in the shadow of the accretion disk around the star 
can continue to accrete, as discussed below. 

In figure S] we show the geometry of the H II region near the point of polar breakout of the ionized gas beyond Vg . In 
this calculation the protostellar evolution has been followed as described in Paper I, but now including the effects of a 
reduction in accretion rate once the H II region breaks out beyond rg. This has only just started to occur at the point 
of the evolution shown in the figure. We have assumed there is negligible reduction in the accretion rate because of the 
Lyman-a radiation feedback since we expect its effects to be limited to relatively small angles around the polar axis. The 
extent of the H II region is calculated using the sector approximation using the density distribution model of Ulrich (1976) 
as described above. We include the effect of electron scattering, but not force due to photoionization, which is discussed 
below. The effect of radiation pressure due to photoionization is strongest for purely radial infall, so its neglect is not 
critical for the models presented here. We do allow for advection of neutrals into the H II region, though they are not very 
important by the time the protostar is ~ 100 Mq. A temperature of 2.5 x 10^ K was adopted for the ionized gas, which 
affects the value of rg. Note that in figure |4] we have assumed an infinitely thin accretion disk. The polar and equatorial 
breakout conditions are shown as a function of /kcp in figure [5] Once the protostar has reached the masses indicated by 
the "Equatorial" line in this figure, we do not expect accretion to be able to proceed from directions that have a direct 
line of sight to the protostar, i.e. those directions that are not shielded by the accretion disk. Thus in most cases we do 
not expect H II region breakout to set the final mass of the star, but rather to cause a decrease in accretion efficiency that 
starts to become important at about 50 Mq in the fiducial case. The actual reduction in accretion efficiency depends on 
the thickness of the accretion disk, to be discussed below (Q. 

We can compare the analytic prediction for H II region breakout (eq. I4ip with our numerical calculation, which for the 
fiducial model (/Kep = 0.5, K' — 1, T4 — 2.5, fd = 1/3) finds breakout in the polar direction at a mass of 45.3 Mq. At 
this point the total H-ionizing luminosity is 549 — 2.78 so that (ps = 1-15 and the bolometric luminosity is 5.95 x 10^ Lq 
so that (pEdd — 0.59. With these values, the analytic estimate for the mass at which polar (/i = 1) H II region breakout 
(fHii = fg) occurs (eq. HT|) is 44.7 Mq, in excellent agreement with the numerical value. In Fig. Owe see that the mass at 
which the H II region breaks out does not scale as a simple power of /kcp- This is because 0Edd and 0s vary with stellar 
mass, especially for m* < lOOM©. 



5.2. Radiation Pressure due to Photoionization 

Continuum radiation is dynamically coupled to the gas in the H II region, both through Thomson scattering and through 
photoionization. Since the H II region is optically thin to Thomson scattering, it effectively reduces the force of gravity 
by a factor 0Edd = 1 — L/I^Edd, which as discussed above is ^ 0.2 — 0.4 for to*^2 ~ 1- At distances large enough that 
the mass acting on the gas is to*^, L/LEdd is reduced by a factor 1 + fd. Keep in mind that the decrease in the effective 
gravity due to electron scattering reduces the accretion rate of ionized gas only; it does not affect the accretion of neutral 
gas outside the H II region onto the disk. 
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Fig. 4. — Geometry of the H II region (shaded) assuming the sector approximation (see text) during the breakout phase for the fiducial 
model with K' = 1 and /kgp = 0-5- The protostar is at (0,0) and the disk is in the ^ = plane. At this stage, the star has m« = 45Mq, a 
total ionizing photon luminosity of 549,101 = 2.78, of which S'4g,acc = 0.50 comes from accretion. The H II region has just recently expanded 
beyond rg (at 94 AU) in the polar direction. The geometry of several accretion streamlines is shown by the dashed lines. 
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Fig. 5. — Mass scales of H II region breakout as a function of the rotation parameter /kcp- The lower dashed line marked "Polar" shows 
the mass scale of the protostar at which the H II region reaches Vg (based on star plus disk mass) along the rotation axis of the protostar. 
The upper dashed line marked "Equatorial" shows the mass scale of the protostar when the H II region reaches Vg in a direction just above 
the disk plane (0.97r/2 from the rotation axis). Note this condition for equatorial H II region breakout ignores the effects of reduced accretion 
rates from prior polar H II region breakout, although such effects arc accounted for in the full feedback+evolution models presented below. 
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Every photoionization results in a transfer of momentum hi>i /c to the gas, where hi>i is the mean energy of the photons 
that ionize the gas. The importance of radiation pressure associated with photoionization has long been appreciated 
in studies of active galactic nuclei and X-ray sources (e.g., Tarter & McKee 1973); Haehnelt (1995) pointed out its 
importance in the formation of the first galaxies, and Omukai & Inutsuka (2002) discussed its role in the formation of 
the first stars. They showed that, under the assumption of perfect spherical symmetry, radiation pressure would have 
the counter-intuitive effect of reducing the size of the H II region, thereby eliminating any feedback effect on the growth 
of the protostar. Since photoionizations are balanced by recombinations, the radiative force is given by a^'^^n'^{hi'i/c). 
Generalizing their treatment to include electron scattering, we find that this radiative force balances the effective gravity 



( ^\ = '^EddPc.Gm..^ ^^^^ 



at a critical density ricr given by 



so that ^ ^ 

nc. = 2.15 X 10%dd (^ll) cm-3, (43) 

where we have assumed a typical ionizing photon energy of 1.5 x 13.6 eV. For gas accreting in free-fall (i.e., it enters the 
H II region at a velocity that is unaffected by radiation pressure) , this corresponds to a critical radius 

r„ = 2.36xlOiVlddh^4 cm, (44) 



2.5 



"Sd, -3 
1-6 ^27/7 



-5.49xlOi34dd(||j cm. (45) 

Even for 0Edd = 1, this radius is typically a few AU in size, and the infall velocity is highly supersonic relative to the 
sound speed of the ionized gas. Omukai & Inutsuka showed that as the ionizing flux from the protostar increased, the 
radius of the H II region would increase until it approached rcr. As it did so, the inflow velocity would drop, the density 
would rise, and the accreting gas would be able to absorb a larger number of ionizing photons. In fact, the accretion flow 
inside could absorb more ionizing photons than any star could emit. As a result, the H II region would remain trapped 
at small radii and the gas would continue to accrete supersonically onto the star. It is not clear that this flow would be 
stable in three dimensions, however, since a fluctuation that placed an element of ionized gas at r > would result in a 
net outward force on the gas. 

Angular momentum in the accreting gas changes this picture completely: The density of the infalling gas is reduced 
inside the disk radius (eq. I16p . which is generally much larger than rcr: 

For typical masses ~ IO^Mq, the disk radius is larger than only for very low rotation, /kcp ^ O.O20Edd- Repeating 
the analysis that led to equation with the density appropriate for rotating infall (eqs. and . we find 

2V>(l±li\'' ryi__^__,r^l\ (47) 

where rcr, spherical is the critical radius for spherical infall (eq. I45p . This result shows that in a rotating infall, the critical 
radius beyond which radiation pressure dominates the effective gravity is intermediate between the critical radius in the 
spherical case and the disk radius. Comparison with equation (|39p shows that in the rotating case, the critical radius 
is comparable to the gravitational radius, where pressure effects can drive the outward expansion of the H II region. It 
therefore appears that for typical values of the rotation, radiation pressure due to photoionization cannot result in the 
confinement of the H II region; correspondingly, feedback by the H II region cannot be curtailed by this effect. 
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6. DISK SHADOWING 

An optically thick accretion disk is able to shield part of the outer accretion flow from direct protostellar feedback. In 
order to determine how effective this shielding is, we must know the thickness of the accretion disk. With a few significant 
exceptions (e.g., Paczynski & Bisnovati-Kogan 1981; Meyer & Meyer-Hofmeister 1982; D'Alessio et al. 1998), almost all 
the work on accretion disks has gone into determining their radial structure; the vertical structure is generally integrated 
over. Here we shall estimate the thickness of the disk under the assumption that it is geometrically thin but very optically 
thick. We neglect convective and turbulent transport in the disk, which D'Alessio et al. found to be small for the cases 
they considered. We focus on the inner parts of the disk, where self-gravity is unimportant (Tan & Blackman 2004) . For 
simplicity, we neglect heating of the disk by irradiation from the central source; our estimate of the disk thickness is thus 
a lower limit to the true thickness. 
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The radial structure of the disk is governed by the equations of energy conservation and of angular momentum conser- 
vation. Energy conservation gives the emergent flux as (Paper I) 



Here 

/n7o\ 1/2 



(50) 



/^l- — , (51) 

is the factor that embodies the boundary condition that angular momentum cannot be transferred across a surface on 
which the angular velocity has no gradient; vjq is the cylindrical radius at which dfl/dw vanishes, which we take to be 
equal to the stellar radius. The dimensionless factor 0/ describes the advection of thermal and internal energy in the disk 
and is generally less than unity. 

To evaluate the angular momentum transfer in the disk, we adopt the a-disk model of Shakura & Sunyaev (1973), in 
which the transverse stress in the disk is proportional to the pressure, — — (we have included the factor | to 
conform with convention — Frank et al. 1995). The equation describing angular momentum transport is then 



m^nf = Gna / Pdz, (52) 
Jo 

where Zg is the height of the surface of the disk. 

The vertical structure of the disk is governed by three equations: First is the first moment of the radiative transfer 
equation, 

~a — = ' (^^) 

oz c 

where Kp is the flux-weighted mean opacity per unit mass and F{z) is the radiative flux. We assume that the effective 
optical depth for true absorption, t* — (TabsTscatt)^^^, is much greater than unity so that the gas is approximately in LTE 
(Shakura & Sunyaev 1973; Artemova et al. 1996). Then P^ad — ^o,T"^ ^nd Kp ~ kr, where is the Rosseland mean 
opacity per unit mass, so that equation l[53|l reduces to the equation of radiative diffusion, 

dT 3krpF 



dz 16aT3 ■ 

The second equation describes the growth of the flux due to viscous dissipation, 

OF , dn 9 



(54) 



— = -(piw^^w —— = -(piaQP (55) 
oz ovj 4 

(Shu 1992). We have included the factor 0/ to allow for the reduction in the flux by the advection of internal energy. 
In addition to the factor |(/)/, equation ([55)1 differs from the expression adopted by Shakura & Sunyaev (1973) in that it 
has dF/dz oc P rather than cx p. One can show, however, that the height of the disk is very insensitive to this change. 
Integration of equation ((55)) together with equation l[52|) leads directly to the energy equation ((50)) . 
Finally, we have the equation of hydrostatic equilibrium, 

dP pGm^z 
oz 



P = Pg+Pr,d = — + ^^- (57) 



where the pressure P includes both gas pressure and radiation pressure, 

pkT AaT^ 

p 6C 

For a primordial gas with a helium fraction of 0.079, the mean mass per particle is 

1.32mH 2.20 x lO'^* 

^"t:o8T^= 1.08+x. 

where Xe = rie/n^ is the ionization fraction relative to hydrogen; for a fully ionized primordial gas, p = 0.98 x 10^^*^ g. 
With the aid of equation (I54p , the equation of hydrostatic equilibrium becomes 



OZ \ p J 



Gm^pz pkrF / 3kcp \ 



(59) 



This is a two-point boundary value problem, in which F = Fq at the surface of the disk and F = at the midplane. We 
assume that the disk is very opaque, so that we can neglect the flux generated above the photosphere; we can therefore 
apply the surface boundary conditions at the photosphere, located at Zph. Since we have assumed that the disk is opaque 
and are neglecting irradiation, the surface temperature is small compared to the central temperature; as a result, the scale 
height near the surface is small and Zph ~ Zs- The temperature at the photosphere is the effective temperature, which is 
related to the emergent flux by Fq = aT^g. 
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We have addressed this problem both analytically and numerically. The case of pure radiation pressure is trivial to 
treat analytically, since hydrostatic equilibrium gives 

KnF Gm^z 

= 5- (60 

c w"^ 

which is true throughout the disk. At the surface of the disk (which is denoted Zsr for a radiation-pressure supported 
disk) this gives 



Gra^,c ' 

where kr_s is the opacity at the disk surface. With the aid of equation (|50l). this becomes 

/ 30/ kr. 



(61) 



V 87 



TO,/ (62) 



8.77 X 10^° </./TO*,_3/ cm, (63) 

where we assume that m, = m*d/ (1 + fd) 1^,^ and where kt — ("■e/p)o'T — 0.35 g cm^^ is the opacity due to electron 
scattering for fully ionized primordial gas. Note that the thickness of a radiation-supported disk depends only weakly on 
radius through the factor / = 1 — (tzjo/w)^/^ and possibly through a variation in the opacity at the surface. 

The case of a disk supported by gas pressure is more complicated, and is discussed in Appendix E for the case of 
constant opacity. There we show that the height of such a disk is 

/ , - X 1/10 / X 21/20 , . j.xl/5 

^..-1.21xl0-(ii^) (f) (I^^:^^ cm, (64) 

where we have normalized a to a typical value of 0.01. In the above expression, we have replaced the constant opacity in 
the derivation in the Appendix with a suitable mean value. Observe that the height of a disk supported by gas pressure 
scales almost linearly with tn, so that the aspect ratio is approximately constant. 

Numerical solution of the structure equations shows that when both radiation pressure and gas pressure are important, 
the approximation 

z.,^(z%^ + z%^)'" (65) 

is accurate to better than 10% over the range 10^'' a < 10^^, kr ~ kt, 10'' K< Toff 10^ K, and 0.01 < Zsr/zsg < 10, 
provided the opacity is constant. 

We have also solved the equations for the vertical structure of the disk numerically with a realistic opacity variation 
with density and temperature (Iglesias & Rogers 1996). We follow the disk structure during the course of the protostellar 
evolution (i.e. as to,, to, and r.^, evolve). We adopt an a-viscosity parameter of 0.01, typical of values associated with the 
magnetorotational instability (Balbus & Hawley 1998; Tan & Blackman 2004), although the disk thickness is not very 
sensitive to this choice. An example of the vertical disk structure is shown in figure [S] for a location at lOi?, around a 
IOOM0 main sequence star, accreting at 2.4 x 10""^ Mq yr"-', i.e. the fiducial rate from a. K' = 1 core with no reduction 
due to feedback. The numerical value of Zs/tu is 0.33. This compares with an analytic estimate based on equation 
of 0.31, where we adopted kr ^ = 0.75 cm^ g^^ and kr = 0.6 cm^ g~^ (Fig- E]). 

The radial variation in the density scaleheight and disk surface height for the above model is shown in figure [T] Note 
that Zs/w is approximately constant with w. For simplicity, in our numerical models we evaluate Za/vj at a radial location 
Tu — lOi?, and use this to evaluate the fraction of the sky shielded by the disk. In the example shown in figures (O and 
(O, Zs/tu — 0.33 at this location, and the fraction of the sky shadowed by this disk photosphere is fsh — 0.31. If the 
matter at infinity were spherically distributed in the envelope, then this would be the approximate accretion efficiency 
once the H II region had expanded to large distances (and assuming all material in the shadow of the disk remained 
neutral). Note that the disk model in the numerical example given above is somewhat thicker than would be present 
around a 100 solar mass star accreting from a K' = 1 core, since the accretion rate would have been reduced by feedback. 
Our numerical models of the growth of the protostar account for such effects self-consistently. 



7. DISK PHOTOEVAPORATION 

As we have seen, in the presence of rotation the various feedback mechanisms discussed above will first disrupt infall 
in the polar direction and may leave behind much of the material in the equatorial plane. Material close to the plane 
will be shielded from the feedback effects by the formation of an accretion disk. However, this gas is subject to disk 
photoevaporation, and accretion will cease when the photoevaporation rate exceeds the accretion rate onto the star-disk 
system. 

To estimate when this criterion is reached, we apply the model of HoUenbach et al. (1994) to estimate the rate of 
photoevaporation from the disk. This rate is calculated assuming a steady disk with no infall from above or below. The 
diffuse ionizing flux, reprocessed through the flared atmosphere of the disk, illuminates and ionizes material near and 
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Fic. 6. Vertical structure of the accretion disk at r = lOr* ~ 43Rq for m* = 100 Mq, K' = 1, /kbp = 0.5 and no reduction in accretion 
efficiency due to feedback. 
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Fig. 7. — Radial dependence of the aspect ratio of the disk surface, Zs/w, and disk density vertical scaleheight, h/vj, for m* = 100 Mq, 
K' = 1, /kbp = 0.5 and no reduction in accretion efficiency due to feedback. 



beyond rg. HoUcnbach ct al. considered the possibility that the disk would be flattened due to a stellar wind, but we 
shall use the results of their weak wind case. The photoevaporation rate is calculated via 

"levap = / 2-n:no{r)r dr, (66) 
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where the flow velocity, v, is set equal to the ionized gas sound speed, q = 18.4(T4/2.5)^/^ km s ^, and uq is the density 
of ionized gas at the base of the ionized disk atmosphere. Their analysis gives 

mcvap = 4.1 X 10-' Sli'^T^-Wjf^^ Mq yr-i. (67) 

As they acknowledge, this result is quite approximate, and a numerical study of this problem would be worthwhile. For 
primordial stars with an ionizing luminosity given by equation ([34]), the photoevaporation rate becomes 

m,vap = 1.70 X 10-^5^^(1 + ml% Mq yr-\ (68) 

There are two corrections that could be applied to this result. First, Begelman, McKee, & Shields (1983) and Woods 
et al. (1996) showed that for analogous winds from AGN disks, the flow can start from radii well inside Vg. Numerically 
integrating the expression given by Woods et al. , we find that mass loss inside rg increases the total mass loss by a factor 
1.5. On the other hand, radiation pressure due to electron scattering reduces the effective mass by a factor (fiEdd ^0.3 
(eq. HD|) . Since the mass loss rate scales as the square root of the gravitational mass, this reduction approximately cancels 
the increase due to mass loss from the inner disk. We therefore adopt the HoUenbach et al. estimate of mcvap for our 
analytic and numerical estimates. 

Accretion onto the star will cease shortly after the photoevaporation rate exceeds the accretion rate onto the star-disk 
system, which is given by equation From equation ([55)1 . we find that the resulting maximum stellar mass is 

28/47 12/47^,60/47 /r,t-\0.24 

Maxm,^^,^6.3 . (69) 



vf-f (l + /<i)26/47 KT,^ 

Recall that e*(j is the instantaneous star formation efficiency — i.e., the ratio of the accretion rate onto the star to the rate 
that would have occurred in the absence of feedback. In the present case, this ratio is just the shadowing factor, /sh, 
introduced in the last section. In the numerical solution described below, we keep track of e*^ as a function of time; for 
the analytic case, we make the simple approximation that the shadowing sets in when the stellar mass reaches mi , so that 
= 1 until the mass of the central star reaches mi and e,^ — fsh thereafter. It is then straightforward to show that 

e*d = -j Tj 7 7, (70) 

provided that m^^ > mi. Note that the average efficiency = 1 at the onset of shadowing (m»rf = mi) and that 
e*d fsh at late times m^^ ^ mi. Normalizing fsh to a typical value of 0.2 from the results of ^and allowing for smaller 
accretion rates due to feedback, we find 

Ma. 1,45 /.»/"(!) (ii) , (71) 

where we have set the ionizing luminosity factor (j)s — I and the disk mass fraction fd = -^'j we have normalized e*d to a 
value of 0.25, which is approximately correct for K' — 1 and for mi ~ 5OM0 as found in fe. 1.21 and m*c( = 200Mq. This 
analytic estimate therefore suggests that for the fiducial case {K' — 1) the mass of a first generation star should be of 
order 140Mq. We now confirm this with more accurate numerical integrations. 

We evaluate the photoevaporative mass loss rate in our numerical model with feedback (Figure [5]) . The accretion rate 
is reduced as the H II region, with T — 2.5 x 10^ K, expands to rg and beyond, although accretion is assumed to continue 
from directions shielded by the photosphere of the accretion disk. The disk and protostellar structure and feedback are 
calculated self-consistently given this evolution in the accretion rate. Beyond about 45 Mq the accretion efficiency starts 
to be reduced below unity. By about 137 Mq the evaporative mass loss rate has become greater than the accretion rate 
and we then expect very limited growth of the protostar. We identify this mass scale as our fiducial estimate for the 
initial mass of the first stars. At this stage fsh — 0.19 and ^5 = 1.370 This estimate of the mass at which accretion end 
agrees well with the analytic estimate of equation ([7T|l . 

We investigate the sensitivity of this result to the ionized gas temperature by setting this equal to 50,000 K (Fig. [9|). 
This causes the H II region to break out sooner and for the disk photoevaporative mass loss rate to be higher. However, 
with the other parameters unchanged {K' = 1, /kop = 0.5), this has only a modest effect on the final mass, reducing it 
from 137 Mq to 120 Mq. 

We also consider the effect of changing the entropy parameter of the initial core by factors of two to higher and lower 
values {K' = 0.5,2) (Fig. [9]). This corresponds to a change in accretion rate of factors of 4.4 since m* cx (X')^^/^. H II 
region breakout is accelerated/delayed by about a factor of two in protostellar mass by these changes. The final stellar 
mass set by disk photoevaporation shows a slightly broader range of factors of 2.4 smaller/greater than the fiducial value. 
This is consistent with equation which would predict a change of a factor of 2.4 if the disk thickness was constant 
and accretion ionizing luminosity negligible. 

Finally we explore the effect of changing the core rotation. FigurefTUlshows models with /xop = 0.0625, 0.125, 0.25, 0.5, 0.75 
for K' = 1 and Ti = 25, 000 K. Cores with higher rotation rates have lower densities in the infall envelope near the star 

Note that only about 10% of this excess H-ionizing photon production rate is due to accretion. The remainder is due to our assumption the 
spectrum of the protostar can be approximated as a blackbody, rather than the more detailed stellar atmosphere models of Schaerer (2002). 
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so the H II region can break out more easily. However, for rotation parameters /kcp ^ 0.25 this makes little difference 
to the final mass, which is set by the balance between (inner) disk-shadowed accretion and photoevaporative mass loss. 
For smaller rotation parameters (/kcp ^ 0.125) the process of H II region breakout does play an important role in setting 
the mass scale at which the accretion rate is truncated to be smaller than the photoevaporative mass loss rate. However, 
given the results of numerical simulations of primordial core formation (O'Shea & Norman 2007), it appears that these 
low values of rotation are very rare. 
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Fig. 8. — Feedback-limited accretion: fiducial case. The evolution of the accretion rate versus protostellar mass is shown for the fiducial 
model (/kop = 0.5, K' = 1, Ti = 25,000 K) in the cases of "no feedback" and "with feedback". In the latter, the accretion efficiency is 
reduced as the H II region expands to Vg and beyond. However, accretion is allowed to continue from directions that are shadowed by the 
disk photosphere. The disk structure and protostellar structure and feedback are calculated self-consistently given the evolution in m». Also 
shown is the photoevaporative mass loss rate, riicvap, which starts once the H II region has broken out in the equatorial direction and grows 
as the ionizing flux increases. We see that this mass loss rate becomes greater than the accretion rate at m* ~ 137 Mq, and we identify this 
mass scale as our best estimate of initial mass scale of the first stars. 



8. CONCLUSIONS 

Recent numerical studies have indicated that the initial conditions for primordial star formation are dense, massive gas 
cores in approximate hydrostatic and virial equilibrium. These physical properties are set mostly by the microphysics of 
H2 cooling and not by the initial cosmological density perturbations. 

We have described the rate of collapse of these gas cores as a function of the entropy parameter of the gas, and the 
amount of mass that already collapsed. This accretion rate is very large, so that once an optically thick protostellar core 
forms the star grows very quickly. 

We have developed a simplified method for modeling protostellar evolution and applied the appropriate accretion rate 
for primordial protostars. The method allows for the treatment of accretion of gas with angular momentum, so that part of 
the accretion occurs via a disk. Using a realistic degree of rotation for the initial gas core we find that, after the protostar 
has grown to about a solar mass, essentially all of the accretion flow is via the disk and conditions at the protostar are 
optically thin, in contrast to the spherical case. This means that the radiation field that the accretion envelope is exposed 
to is significantly hotter so that ionization and FUV radiation feedback can become important. 

We considered the impact of the protostellar feedback on the infalling envelope. Again rotation is important because 
it modifies the density distribution in the vicinity of the star. First we discussed the effects of photodissociation of H2, 
the primary coolant. We showed that this does not stop accretion if the protostar has already begun to form, but it 
can suppress star formation in the vicinity (c.f. Glover & Brand 2001). Next we considered radiation pressure feedback 
due to resonant scattering of FUV radiation in the Lyman-a damping wings. As a result of the high column densities of 
neutral gas around the H II region, this radiation is trapped and the pressure amplified by large factors. This radiation 
pressure becomes larger than the ram pressure of the infalling gas in the polar directions for stellar masses of order 2OM0. 
However, once the infall is reversed at the poles, the Lyman-a photons can escape and the accretion in other directions 
proceeds unimpeded. We then considered the growth of the ionized region. Once the expansion velocity of this region 
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Fig. 9. — Feedback-limited accretion: effect of ionized gas temperature and accretion rate. The fiducial model (/kcp = 0.5, K' = 1, 
Ti_4 = 2.5 K) shown in Fig.[8]is compared to models in which one parameter has been changed: a model with Ti^4 = 5 and two models with 
K' = 0.5,2. The dashed lines show the accretion rate to the star, rri*, and the solid lines show the photoevaporative mass loss rate, rhevap. 
The change in temperature causes relatively minor differences, while the change in K', equivalent to a change in m« of factors of 4.4 above 
and below the fiducial level, leads to roughly a factor of 2.4 change in the final stellar mass. Note the increase in m* for the K' = 0.5 case 
at around 35Mq is due to a thickening of the inner accretion disk as the star contracts down to its main sequence configuration and assumes 
material at large distances still remains to be accreted in the enlarged shadowed region. 



exceeds the free-fall velocity, the accretion is halted. This typically occurred at about 50 to lOOM©, although it took much 
larger masses for cases with little angular momentum. The ionized gas is confined to the region above the disk, however, 
so accretion can continue in the shadow of the disk. Evaluating this, we found that shadowing permitted accretion to 
continue at a rate of about 20-30% of that in the absence of the H II region. Allowing for photoevaporation of the disk, 
we found that the final stellar mass is about 140Mq in the fiducial case. 

Table [T] summarizes how the mass scales set by protostellar feedback depend on model parameters. The final mass of 
a Pop III.l star depends fairly sensitively on the entropy parameter of the accreting gas (i.e., approximately as {K')^-^), 
which in turn determines the overall accretion rate to the star-|-disk, but not very much on core rotation (for /Kep 0.25) 
or ionized gas temperature (Ti). At very low values of core rotation, H II region breakout is delayed until high protostellar 
masses, at which point the disk photoevaporation rate soon exceeds the residual disk-shadowed accretion rate. However, 
these small values of /xep are not very likely to occur in nature. 

The final masses predicted by our model overlap the range of masses necessary to produce pair instability supernovae, 
140 — 260 Mq (Heger & Woosley 2002). Rotation may lower these limits (S. Woosley, private comm.). The lack of the 
expected nucleosynthetic signature of such supernovae in the abundance patterns of very metal poor halo stars (Tumlinson 
et al. 2004), could indicate that such massive Pop III.l stars were relatively rare or that they tended to enrich regions not 
probed by typical halo stars, perhaps the centers of larger galactic halos. The conclusion by Scannapieco et al. (2006) 
that Pop III star formation should be fairly widespread in regions now probed by Galactic halo stars, can be reconciled 
with the abundance pattern observations if most of this star formation leads to either Pop III.l stars from relatively low 
entropy {K' < 1) gas cores or Pop III. 2 stars that also have a mass scale below the pair instability threshold (see also the 
study by Greif & Bromm 2006). Further work is required to determine the range of pre-stellar core parameters, primarily 
K' and /kgp, exhibited in cosmological simulations, in order to predict the frequency of pair instability supernovae. 

One may ask how the feedback mechanisms we have considered relate to those that operate in contemporary massive 
star formation. We note that the maximum mass attained in our fiducial model of PopIII.l star formation is very similar 
to that inferred observationally in local massive star clusters (e.g. Figer 2005). However, after decades of study, it remains 
unclear whether the maximum mass of stars forming today is set by feedback or instabilities in very massive stars (Larson 
& Starrfield 1971). Here, we have argued that the maximum mass of primordial stars is set by feedback. The primary 
differences in the feedback processes then and now are: 

(1) Dust. In contemporary star-forming regions, dust destroys Lyman-a photons, eliminating them as a significant 
pressure. On the other hand, the dust couples the pressure of the UV continuum radiation to the gas very effectively. 
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Fig. 10. — Feedback-limited accretion: effect of rotation. The fiducial model (/kcp = 0.5, K' = 1, Ti.4 = 2.5 K) shown in Fig.[8]is compared 
to models in which only the rotation parameter /kcp has been changed: /kcp = 0.0625,0.125,0.25,0.75. Smaller rotation parameters result 
in higher polar gas densities in the infall envelope and thus delayed H II region breakout (Fig. [Sjl. However for /kcp 0.25 this has relatively 
little effect on the final mass, which is set by disk photoevaporation (note the convergence of the /Kep = 0.25,0.5,0.75 models). At smaller 
rotation parameters the process of H II region breakout plays an important role in setting the mass scale at which the accretion rate is 
truncated to be smaller than the photoevaporative mass loss rate. 



and it remains to be determined whether this hmits the final mass of the star; e.g., Yorke & Sonnhalter (2002) find that 
it does, whereas Krumholz et al. (2005a) have not found evidence that it does. Dust also affects the evolution of H II 
regions, absorbing a significant fraction of the ionizing photons in dense H II regions (Spitzer 1978), thereby reducing the 



Table 1 

Mass Scales of Population III.l Protostellar Feedback 



K' 


/kcp 


T,/(104 K) 






nT-*,ovap {MqY 


1 
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^Mass scale of H II region polar breakout. 

'^Mass scale of H II region near-equatorial breakout. 

'^Mass scale of disk photoevaporation limited accretion. 

'^Fiducial model. 

®This mass is greater than mt^evap in this case because it is calculated without allowing 
for a reduction in m, during the evolution due to polar H II region breakout. 
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H II feedback discussed in §5 and the photoevaporation in §7. 

(2) Magnetic fields. In contemporary protostars, magnetic fields drive powerful protostellar winds that drive away a 
significant fraction of the core out of which the star is forming (Matzner & McKee 2000). The cavities created by these 
winds allow radiation to escape from the vicinity of the protostar, significantly reducing the radiation pressure (Krumholz 
ct al. 2GG5b). Tan & Blackman (2004) considered the influence of such outflows on Pop III.l cores, concluding the 
instantaneous efficiency of accretion could be reduced by a factor of about 2 from the no-feedback case in an isotropic core 
by the time the star reached lOOM©. However, these outflows would not confine ionizing feedback from the star at these 
masses, so much of the gas that could be expelled by outflows would have already been disrupted by H II region breakout. 
We conclude that outflows would have a relatively minor effect on the results presented here, and that the masses of the 
first stars arc mostly influenced by radiative fcc;clbac;k. See Tan & McKee (2008) for further discission. 

(3) Stellar temperatures and luminosities. Primordial stars were significantly hotter than contemporary stars, resulting 
in significantly greater ionizing luminosities. In addition, the accretion rates of primordial massive stars are much greater, 
at least initially, than those of contemporary massive stars (McKec & Tan 2002; 2003). Future calculations will show 
whether feedback can be as effective in setting the maximum mass of contemporary stars as we have argued that it is for 
primordial stars. 
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APPENDIX 

A. LINE PROFILE WITH RAYLEIGH SCATTERING 

A complication that occurs in our problem is that the column density can become so largo that Rayleigh scattering is 
important. Prom Jackson (1975), we find that in general the scattering cross section can be expressed as 

a{i') = aAuDCpiu (Al) 

where 

^ = / "^^""^ ^^^^ 

n (4^V^g)(r/47 



and r is the total spontaneous transition rate out of the upper and lower levels of the transition. For the simple case of a 

two-level atom (which can be used for Lyman a), T — A21 is the Einstein A coefficient for the transition. Physically, (pi^, 
is the line profile for scattering by an individual atom; the line profile for a gas, (pu, is obtained by convolving <j)ii, with a 
Maxwcllian distribution. 

We now develop an accurate approximation for Defining 

4i/4 



i\ im^).m ^^^^ 



n (r/4.)/(.) (^g^ 



which is unity at line center, we have 



+ (r/47r)2' 

The approximation of dropping the factor {y / Vq)'^ f {y) in the denominator has an error of order r/(47r!/o) = A2i/(47rfo), 
which is less than 10~^ for Lyman a. At low frequencies {v <C J^o), we have cpiu oc f{v) oc v'^, the standard frequency 
scaling for Rayleigh scattering. 

In terms of the normalized frequency shift x = IS.vIIS.vy), wc have cr(x) = a(^\x where 



TT 



cfl -V a;2 



(A7) 



and a = r/(47rA:/£i). For Lyman a, a = 6.04 x 10 ^/Auc6i where Aw^e = Aw£)/(10^ cm s The line profile for a 
gas, (px-, is the same as this in the line wings (x 1). We conclude that the line profile in the damping wings is given by 

(Px^^^ (x»l). (AS) 

The correction to the usual expression for the damping profile is given by the factor f{v), which drops to 0.5 a.t v = O.SOvq. 
The relation between v and x is given by 

— = 1 + ^^-^ = 1 -I- 3.33 X 10-^Avd 6X, (A9) 
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where the final expression is for Lyman a. 

The optical depth at a frequency labeled by x is Ta; = fcs_^x- For Lyman a, 

feff = L34 X 10-"iVcff/AwD,6, (AlO) 

(Neufeld 1990), where iVoff is the effective column density of H I (cf. eq. [^0]) . Thus, in the damping wings of Lyman a we 
have 



T.r = 2580 



(x> 1). 



(All) 



The frequency O.Svq, where /{v) = 0.5, corresponds to a; = GOOO/Auc^g- In order to have optical depth unity at this 
point, the column density must be A'cff ~ 3 x lO'^^. 

B. ENHANCEMENT OF LVMAN-a INTENSITY IN AN OPTICALLY THICK MEDIUM 

Here we derive the increase in Lyman-a intensity relative to the optically thin limit due to the trapping of photons 
in regions of high column densities, such as the neutral gas around the protostellar H II region. This factor was used in 
for the calculation of the Lyman-a radiation pressure feedback, and some of the symbols in this appendix are defined 
there. We first consider the case of pure scattering, and then include the effect of destruction by two-photon emission. 



B.l. Case of Pure Scattering 

We follow the treatment of Neufeld (1990), which extended earlier work by Harrington (1973). He considered a uniform, 
planar slab of thickness 2fi, with the origin at the center of the slab. A planar source of photons located at fg produces 
an incident flux Fi in each direction; he normalized to Fi = 1/2. We assume that there is no absorption. First consider 
the case in which the photons are injected at line center (xi = 0). The frequency-integrated intensity at a point f in the 
slab is 



J(t; Xi = 0) 
2F, 



1/2 



4aTr 



1/3 r 



2fL 



2fL 



2fr 



where 



In general. 



-y6r(i) cos(n7ru)) 



127r7/3 

n— 1 



1 



,4/3 



3TL0i 



(Bl) 



(B2) 



(B3) 



where — 4>{xi) is the line profile at the injection frequency; we have assumed that the photons are injected at line 
center, so that (pi = 0o is not small. We shall further assume that — Ts, which is valid provided A-fs = — fg 3> 1 — i.e., 
the source is not too near the edge of the slab. 

For the case in which the source is at the center of the slab {fg = 0), the intensity at the center of the slab is 

"r(i)C(|)- 



Js{x,=Q) = J{Q: = 0)= 24/3_i 



{afLf'^'F, = 0.4362(afi)i/3F„ 



(B4) 



22/3^8/3 

as originally found by Harrington (1973). This has the same scaling as expected from the heuristic argument given in 
m in the absence of any scattering, the mean intensity would be J = Fi/27r (assuming isotropic emission), so the mean 
intensity is indeed enhanced by a factor of order (afi)^/^. 

Next, consider the case in which the source is near the edge of the slab, but, in view of the smallness of the Lyman-a 
mean free path, still at a large optical depth from the edge (f 3> Af^ 3> 1). The maximum intensity occurs at -fg, and is 
proportional to 

'AfA V6r(l) - 1 " cos (li^ 
TL ) 127rV3 2^ 



74/3 



(B5) 



Approximating the sum by 



dn n "'/^ sin^ 



rAi 



2fL 



COS,- 



IT At 

TL 



1/3 



(B6) 



(Gradshteyn & Ryzhik 1965), we find 

Js{x, = 0) = 0.518(aAr,)i/^F, = OAn{af,B)^/^ F,, (B7) 

where the effective optical depth is Toff = 2A-fs (eq. [20)1 for the case in which the source is close to one edge. Note that 
this agrees quite well with our ansatz, equation l|20p . since when the result for a source near the edge is expressed in terms 
of feff, it is nearly the same as that for a source near the center (eq. IB4p . 
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Neufeld's results are valid only in the limit in which utl ^ 10^, so that the transfer is completely determined by the 
damping wings. Hummer & Kunasz (1980) show that for smaller optical depths, there is an intermediate regime in which 
J/Fi is about constant. They define a quantity pjjK = PkD/tl, where the density p is assumed to be constant and where 
D is the mean distance traveled by an escaping photon. As shown by Ivanov (1970), this is (47r/2fiFi) J Jdf, which 
in turn is proportional to Js/Fi. For the case 1 > Auu^g ^ 0.1 (corresponding to 10^ K> T > 10^ K), the results of 
Hummer & Kunasz (1980) imply that phk, and therefore Js/Fi, are within about 0.1 dex of their values at af^ = 450 for 
10^ > qtl ^ 1. We assume that the same behavior obtains in terms of foff for a source near the edge of the slab. Thus, 
for foff ^ 1/a = 1660Aw£)_6i we have 

Js 



F, 



0.411(afeff)V3 

Min [l, (afcff/450)i/3] ' ^ ' 

For 450 > afcs > 1, this gives Js/Fi ~ 3.15. The condition afcs > 1 for the validity of equation (|B8|1 corresponds to to 
column densities A^off ^ 10^^A?;|, g cm^^. 

This result is valid for Lyman-a photons produced by the H II region, since such photons are very near line center. 
Stellar FUV photons are not restricted to line center, but fortunately Ncufcld has evaluated the intensity for an arbitrary 
injection frequency. Since his results for this case are somewhat complicated, so we shall evaluate the intensity far from 
line center and then smoothly join the result onto the result we have found above. If the injection frequency is large 
compared to the diffusion frequency [xi 3> (afoft)^/'^], then the photons scatter approximately coherently. The intensity in 
this case can be found either from Neufeld's general resultf0 or from a simple solution to the radiative transfer equation. 
We normalize the injection frequency, 

Xi = -—^ vlTiT' (^^) 
(aroff)^/'^ 

so that photons injected with ii <C 1 are in the diffusion regime and those with 3> 1 are in the coherent scattering 
regime. In the latter case, we find 

^..iso 3 f 3 ^ (afeff)i/V(z.) 

4^ " 1^4^ ) if ^ ' ^ ^' ^ ^ 

where /(z/) is the Rayleigh-scattering factor defined in equation (|A4|) . We have put the subscript "iso" on the mean 
intensity to indicate that it has been derived under the assumption that the medium is optically thick at the frequency 
Xi so that the photons are isotropic. Sufficiently far in the line wings, equation (jBlOp shows that Js, iso goes to zero. This 
approximation is developed further in Appendix [C] below. If several lines contribute to the opacity at a given frequency, 
then the right hand side of equation (jBlOp should be summed over the lines, since it is the total opacity that governs the 
mean intensity. 

We are now in a position to join the result for the frequency diffusion (eq. IB7[) to that for coherent scattering in the far 
wings of the line (cq. IBlOp . Taking the harmonic mean of these results, we obtain an expression that is approximately 
valid for all injection frequencies: 

J., iso 0.41 l(afeff) 

(Bll) 



F, 



F, 



Min[l, (afeff /450) 1/3] +5.41 [f|//(z.)]- 



Note that the numerical coefficient has been chosen to agree with the case of a source near the edge of a cloud, which is 
the one most relevant to our problem. The intensity is half that at line center for an injection frequency Xi — 0.43 (for 
afoff > 450 and /[z/] ~ 1). 

B.2. Effect of Two-Photon Emission 

In the absence of dust or molecular hydrogen, the dominant destruction processes for Lyman-a photons are coUisional 
de-excitation (which we shall ignore) and two-photon emission following a coUisional transition from the 2p to the 2s 
state. Lyman-a photons can also be destroyed by photoionization out of the n — 2 state, but since another Lyman-a 
photon is created after the ion recombines, the net destruction by this process vanishes. 

Consider a ID slab of gas with a central source of Lyman-a photons. Let e be the probability of two-photon emission 
per scattering. Then in the limit of large optical depth, the mean intensity at the source is (Harrington 1973) 

Js{x^ = 0) = 0.396(a/e)i/3Fi, (B12) 

which is quite close to the result with no absorption (eq. IB4p with the replacement — > 1/e. Thus, two-photon emission 
prevents the mean intensity from increasing once ef^ > 1. We join this result onto the expression for the case in which 
the source is not at the center and there is no absorption (cq. IB8p by writing 

Jsix.^0)^ ";'^.^^°'°rr . (B13) 

' Min[l, (afeff)i/3] +r' ^ ^ 

with 

r = 1.04(efeff)i/3. (B14) 



4 



The numerical coefficient in Neufeld's equation (2.30) is too small by a factor 3, as confirmed by the author (private communication). 
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To determine the destruction probability e, we need to know the population of the 2s state. In statistical equilibrium, 
this is 

n2p 3 \1 + ne,cr/ne/ 

where 

ne,cr = = 1.55 X 10"* cm^^' (B16) 

<l2s2p 

is the critical density for the 2s —^ 2p transition, A2sis — 8.23 s^"'^ is the two-photon emission rate from the 2s state, and 
q2s2p = 5.31 X 10^'* cm'^ s^^ is the collisional rate coefficient for electron and proton collisions from the 2s state to the 
2p state (Osterbrock 1989). Collisional de-excitation from 2s to Is is much slower, and may be neglected in determining 
this population ratio. The probability per scattering of destroying a Lyman-a photon by two-photon emission is then 

n2sA-2sis 4.4 X 10"^ 

n2pA2pls l+ne,cr/"e 

Collisional de-excitation to the Is level competes with this process for 10^ cm~^, but we shall assume that is less 
than this so that we can ignore this process. Inserting this result into equation (jB14p . we find that the factor that gives 
the effect of two-photon emission is 

^ 0.405(iVeff,2o/Avz5,6)'/^ 



(l+ne,cr/ne)^/3 



(B18) 



We see that at high electron densities (rie 3> Ue^cr), two-photon emission reduces the Lyman-a intensity by a factor that 
depends only on the column density and the velocity dispersion. Although the large column densities of H I needed to 
make this process important occur in regions of neutral hydrogen, one can show that photoionization out of the n — 2 
state is generally sufficiently effective that Ue > rie, cr in the H I just outside the H II regions of massive primordial stars. 
As a result, we have 

r~ 0.405(iVcff,2o/Az;i5,6)'/^ (B19) 

This destruction process operates only for photons that can diffuse to the center of the line, which is where most of the 
scatterings take place. For stellar FUV photons, we assume that this occurs only for photons within a frequency range 
(afoff )^/^Ai/£), or \xi\ < 1, so that T depends on xf. 

m) = {l (B20) 

B.3. The Blackhody Constraint 

We have one last constraint to impose: the intensity we calculate must be less than the appropriate blackbody intensity 
B^. The Lyman-a photons produced by the H II region have a complicated line profile that is concentrated in a frequency 
range ~ 2{afg.sY^^ ^vd- For these photons, we require that the intensity in this frequency range be less than that of a 
blackbody at the temperature of the H II region, 

7 A/r- - ^l/3A n irr \ 0.411(afeff)i/3^^ ^ 

J.,H„ = Mm|2(arc,) A..i?.„ (Th„), ^.^^^^ (^^450) V3] F | ' ^^'^^ 
For stellar photons, the intensity is limited by the blackbody intensity at the stellar surface, B^{T^) = F^^,/tt: 

J ^ ^^ ^ Mm[l, (arcff/450)i/3] +5.41[a;f//(zy)] +F(a;j) J 

The subscript "iso" indicates that Ja^.iso is that part of the stellar radiation that has been isotropized by scattering (see 
Appendix [C] below) . The factor in the denominator ensures that although the integral is taken over the entire spectrum 
of the star, it is only that part near the resonance line that contributes. 

C. ESTIMATE OF THE RADIATION PRESSURE 

In opaque media, the force exerted by radiation can be treated as an isotropic pressure, just like gas pressure. If the 
medium is not opaque, however, the radiation pressure tensor is not diagonal, and radiation pressure does not behave like 
gas pressure. For example, in the optically thin limit, the radiation pressure declines as yet the associated force per 
unit volume is negligible. Here we introduce an approximation to the radiation pressure tensor that separates out the 
isotropic part and show that the radiative force is the gradient of this pressure. 

The radiative force per unit volume is 

f,ad - -V • Prad - - P^F, (CI) 
C 

where Prad is the radiation pressure tensor (e.g., Shu 1991); for a frequency-dependent opacity, kF = J n^Y^. At large 
optical depths, the radiation field is nearly isotropic and the radiation pressure tensor becomes diagonal, Prad -Prad.isol- 
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At small optical depths, the radiation is beamed and Prad — > {^/c)FF, where F is a unit vector. (Note that in one 
dimension, the radiation is not purely beamed at small optical depth, so there is a numerical coefficient < 1 in front of F 
in this expression.) We then have for all optical depths 

Prad ^ Prad, isoI+-FF, (C2) 
C 

and the radiative force is 

frad = -V • Prad ^ - VFrad, iso " -FV • F. (C3) 

c 

Here we have omitted the term F • VF, which vanishes in the optically thin limit since then the flux does not change 
direction along a ray, and which is negligible in the optically thick limit even if the flux does change direction. In the case 
of pure scattering, or in a stellar atmosphere, we have V • F = 0, so the radiative force is then 

frad ^ - VPrad, iso - - puF (for V • F 0). (C4) 
c 

This is the case relevant to our problem. On the other hand, in the case of strong absorption, so that F decreases in a 
distance small compared to the radius, one cannot neglect the term FV • F in evaluating the force. 

If F is known, then Prad. iso can be determined from equation (jC4p once a boundary condition is specified. We assume 
that there is a surface to the gas distribution. It follows that fiad cx; k vanishes outside the surface, so that Prad, iso must 
be constant there. Since the radiation pressure vanishes at infinity, it follows that the constant must be zero. We conclude 
that Prad, iso = at the surface of the gas distribution. 

At large optical depths, the energy density of the radiation is Urad — 4,ttJ/c — 3Prad- We therefore define Urad,iso = 
47rJiso/c = 3Prad,iso- Siuce Mrad ^ F/c at small optical depths, we have 

^rad — ^rad.iso 

+ F/C. 

Under the assumptions of spherical symmetry and no net absorption (V • F = 0), our formulation is equivalent to the 
closure approximation recommended by Shu (1991, pp. 43-44). He shows that the equation for the rr component of the 
radiation pressure tensor is 

^^Prad,rr 1 ^ s 1 



dr 



- (3Prad, rr " "rad) = pnF. (C5) 



Shu points out that 3Prad, rr ~ i*rad — 2P/c, sincc at high optical depths F is negligible and Prad = ^^rad/S, whereas 
at low optical depths Prad,rr — Urad — F / c. With this approximation, the left-hand side of the equation becomes 
9Prad,rr/9r + 2F/cr. Our formulation gives Prad,rr = Prad, iso + F/c. Since we have assumed V • F = 0, we have 
dF/dr = —2F/r in the spherically symmetric case, so that equation (jCSp reduces to equation (|C4p . 

Recall that the Lyman-a radiation pressure, Prad, iso in equation (j22p . was derived for a slab geometry. What is the 
appropriate value of the column density A^off to use in this expression in a more realistic geometry? Here we shall 
determine the relation between a slab geometry and a spherical one; in the next Appendix, we generalize this to non- 
spherical geometries. However, it must be borne in mind that the actual geometry of the infalling gas is far more 
complicated than can be represented by a simple analytic model. 

What is the relation between the optical depth in one dimension and that in three dimensions that results in the same 
radiation pressure? Let F cx r"*"'^, with kp = for slab geometry and kp = 2 for spherical geometry. We assume that 
the density in the spherical case can be described by a power law, p cx r~^i' . For a constant opacity per unit mass, we 
therefore have 

9Prad.iso PokPo ( roX^F+^^P 

' ' , (C6) 



dr c 

where tq is a fiducial radius and po = p(''o)j etc. The optical depth from r to infinity in the spherical case is 



for fcp > 1, so that the solution of equation (jC6P is 

Prad,iso=fTr^^V^)— ■ (^8) 

In order for the radiation pressure in the slab to be the same as that in a sphere for the same value of the flux, we require 

For a free-fall density variation, valid for r > r^^ we have kp = 3/2 so that Acff — (2/5)nr; inside r^, kp — 1/2 is a more 
accurate description, so that A^cff — (2/3)nr. 



Since tid = PhkN^s, we conclude that 
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D. ANISOTROPIC OPTICAL DEPTH AND SUPER-CRITICAL ACCRETION 

In order to estimate how Lyman-a photons escape from the H II region around the protostar, we consider the following 
idealized problem: We assume that the radiation fills a cavity bounded by a thin, opaque shell of variable optical depth, 
T(r). In this case, the flux at the surface of the shell is approximately normal to the surface, F ~ n, and the radiation 
energy density is about constant in the interior of the shell. This model will be approximately valid for Lyman-a radiation 
once the H II region separates from the star, since the H II region provides a cavity in which the optical depth due to 
resonance line scattering is relatively small, so that the radiation becomes approximately uniform there. Equation (jC4[) 
then gives 

F/ \ ^-^rad, iso ^ ^-^rad. iso ^ /t^-i \ 

(r) ^ -c ,\ n ~ -— n. (Dl) 

dT[r) T(rj 



Integration over the surface of the shell gives the luminosity 

L = J F - hdA^ cPrad, 
where A is the total area of the shell and feg is the the harmonic mean optical depth: 

For a spherical shell, this simplifies to 



T(r) Teff 



(D3) 



1 




r dA 






r(r)- 


1 




r dn 




" 4^ j 


' r(f) 



(D4) 

We can generalize this treatment to allow for the possibility that the optical depth is small in some directions. Consider 
the extreme case in which t = over a small area SA — i.e., there is a small hole in the shell. The flux emerging from this 
area is tt times the specific intensity, which is the same as the mean intensity J in the cavity. We therefore find 

F{r) =ttJ= = ~ (r = 0), (D5) 

where the last step follows since we have assumed that the average optical depth is large enough that Pjad ^ F/c so that 
Piad — ^rad, iso (uotc that Prad, iso drops near the hole, but that does not affect the average value of Prad, iso since the hole 
is small) . Combining this result with equation (jPip . we write 



cPrad, 



P(r) ^ (D6) 



as an expression that is approximately valid for all t. With 

1 1 r dA 



TeS A J r(r) + I 

equation (|D2p is valid even if the optical depth is small in some directions. As an example, assume that r = over an 
area SA and t — tq ^ 1 elsewhere. Then we have 



^ 6A\ SA 



A J A 

We require foff I in order for our treatment to be valid, and this will be true if both tq ^ 1 and SA/A <i; 1. 
Equations (jPip and (jD2p imply that the flux at any point on the shell is then 



T(r) \A^ 

In our problem, r < f^s near the poles, so the flux there is enhanced over L/A since radiation originally directed at regions 
of high optical depth tends to escape in regions of low optical depth. This is a quantitative expression for the flashlight 
effect (Yorke & Bodenheimer 1999). 

When radiation pressure is acting against gravity, it is convenient to define the critical flux as the fiux that just 
counterbalances gravity, 

f^cit = (DIO) 

where we have assumed spherical symmetry and where k — a/ fj, is the opacity per unit mass and /x is the mean mass per 
particle. The critical luminosity is then 

Lcrit = 47rr^Pcrit = —■ (Dll) 

cr 

If the opacity is due to electron scattering, the critical luminosity is the Eddington limit. We conclude that supercritical 
accretion — i.e., accretion when L > Lcrit — can occur in directions with r > f^s since it is possible for F to be less than 
Pcrit in those directions: 

F{r) feff 



r(r) 



(D12) 
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For example, an accretion disk can produce supercritical accretion since the optical depth in the plane of the disk is much 
larger than that in other directions. 

This argument works well in our problem because the Lyman-a opacity in the central regions is small due to photoion- 
ization, thereby rendering the radiation approximately uniform there. It is more difhcult to create a uniform radiation 
field in the case of electron scattering in an ionized gas, since the opacity per unit mass is constant and it is difficult to 
create a thin, opaque shell around a star. Nonetheless, the effective optical depth in this case is likely to be of order the 
harmonic mean optical depth, just as we have found for our idealized problem. 



E. VERTICAL STRUCTURE OF AN ACCRETION DISK SUPPORTED BY GAS PRESSURE, WITH CONSTANT OPACITY 

Here we determine the height of an accretion disk supported by gas pressure under the assumption that the opacity per 
unit mass, k, is constant. For a disk supported by gas pressure, the equations describing the radiation field in the disk 
(see can be written as: 



(radiative diffusion), where 



is the surface density above a height z; and 



S = 



dF 

dz 



4(7 

pdz 
FoP 



(El) 
(E2) 
(E3) 

(E4) 



(flux generation). Since gas pressure dominates, we have P = pkT / ^, so that this becomes 

dF _ FqT 

" ~Sc(r)' 

where Sc — Jq" pdz is half the total surface density of the disk and (T) is the mass-weighted average temperature in the 
disk. 

To obtain an approximate solution for these two equations, we set 

1/4 



T 



1 + e(E/S,)^ 
1 + e 



(E5) 



where Tc is the central temperature. The parameters e and £ are to be determined; in particular, e is assumed to be small, 
so that 



1 + 4£(S/I], 
1 + 4e 



Inserting this into equation (jEip . we find 



F = 



AaTl 
3t^ 



l + 4(^+l)e(S/S,)^ 



1 + 4e 



(E6) 



(E7) 



where Tc = krSc is the optical depth from the midplane to the surface. Since = at the midplane, where E = Ec, we 
find 

e = -^^. (E8) 



At the surface (E = 0), we have F = Fq, so that 



Fn 



1 



4crT4 



(l + 4e) 3t, 



and 



F = Fo 



1 



E 

e; 



(E9) 



(ElO) 



Inserting this into equation (|E4p and keeping only the leading term in equation (lESp implies £ — 5/4, so that e = —1/9. 

■4 

eft' 



Since Fq = fTrfr, equation (jE9p implies 



12 ' 



1/4 



Tcs. 



(Ell) 



This approach gives a mass- weighted mean temperature (T) — ^Tc from equations (jE4[) and (|E10[) : on the other hand, 
direct integration of equation (|E5|) gives (T) = ^^c. The 6% difference between these estimates for (T) is a measure of 
the accuracy of our approximations. 
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The equation of hydrostatic equihbrium is 

§ = (E12) 
Approximating T ~ rc(S/Ec)^^'', which is typically accurate to better than 10%, we then find 

where the second step follows from equation (jElip . In terms of the Keplerian velocity vk = {go^Y^'^, the scale height is 
hgc/ru = {cgc/vxY. Since p = —<E/dz, equation (|E13p yields the following equation for E: 



Define the characteristic scale height as 



zdY.'. (E15) 



To obtain an approximate solution to this equation, we adopt the following ansatz for S: 



An approximate evaluation of the integral J zdTi' gives 

3 3z 



zdS' ~ z.E ( i + 1^ ) . (E17) 



Integration of equation (|E15|) then gives 



Sc/ Qwhgc \ Zs J \ 2z. 

which is consistent with the ansatz provided that the height of the disk is 

Z,g = i6Vjhg,f^ , (E19) 

where the subscript g indicates that the height is evaluated for the case in which gas pressure dominates. Shakura & 
Sunyaev (1973) show that the height of the disk ~ {cc/vk)^, where Cc is the central isothermal sound speed. Equation 
(|E19|) implies that in a gas-pressure dominated disk. 



In terms of the sound speed at the photosphere, Cg^off — (fc?cff//i)^/^, this is 



(E20) 



1^ = (540r,)^/« (E21) 
w Vk 

Numerical integration of the structure equations shows that the actual height of a gas-pressure dominated disk ranges 
from IMzgg for Tc = 10'* to l.lOz^g for Tc = 10^, so the approximations made in our analytic estimate are reasonably 
good. 

Paczynski & Bisnovati-Kogan (1981) obtained equation (jE20[) through a different argument: They assumed that the 
disk is polytropic, with P ex p^+^Z", and found that 



Zsg _ [2(n+l)]l/2Cgc 



(E22) 



r VK 

They argued that n is likely to be between 1.5 and 3, so that 2{n + 1) is between 5 and 8; they chose 6 as a typical value. 
We emphasize that our derivation depends only on the assumption that the opacity is constant, and is not based on the 
assumption that the gas is polytropic. 

To complete the determination of the height of a gas-pressure supported disk, we must estimate the optical depth 
through half the disk, Tc — kr'^c- Observe that 

where is the central value of the mean molecular weight. According to the discussion below equation (jElip . the average 
temperature is (T) ~ |Tc, and we adopt this value here. The fact that (T) is so close to Tc justifies setting the mean 
molecular weight equal to /ic, as we have done. From equation (|52p . we then find 

T(. — T . 

QTrai^kTc/ He) 



(E24) 
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Equation (jElip then implies 



2k 

Tcff 



1/5 



_2887ra(fcTeff///c). 

Using this resuh in equation (jE20p for the height of a gas-pressure supported disk, we find 



= 1.21 X 10' 



where we have normahzed a to a typical value of 0.01 



1/10 



21/20 



7/20 

"^*;2 



cm, 



(E25) 



(E26) 
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